Commit e0001b5b authored by Nathan Coppin's avatar Nathan Coppin
Browse files

exploding solved

parent ff649b4c
Pipeline #6488 passed with stage
in 2 minutes and 12 seconds
......@@ -5,7 +5,10 @@ def _advance_particles(particles, f, dt, min_nsub,contact_tol) :
f = np.zeros_like(particles.velocity())
# Compute free velocities of the grains
vn = particles.velocity() + f*dt / particles.mass()
vmax = np.max(np.hypot(vn[:, 0], vn[:, 1]))
if particles.velocity().shape[1] == 2:
vmax = np.max(np.hypot(vn[:, 0], vn[:, 1]))
else:
vmax = np.max(np.sqrt(vn[:, 0]**2+ vn[:, 1]**2 + vn[:,2]**2))
# Estimation of the solid sub-time steps
nsub = max(min_nsub, int(np.ceil((vmax * dt * 8)/min(particles.r()))))
print("NSUB", nsub,"VMAX",vmax, "VMAX * dt", vmax * dt, "r", min(particles.r()))
......
......@@ -25,6 +25,10 @@ set(SRC
include_directories(. ../tools)
add_library(scontact3f SHARED ${SRC})
target_compile_definitions(scontact3f PUBLIC -DDIMENSION=3 -DFRICTION_ENABLED=1 -DROTATION_ENABLED=1)
target_link_libraries(scontact3f mbtools3)
add_library(scontact2 SHARED ${SRC})
target_compile_definitions(scontact2 PUBLIC "-DDIMENSION=2")
target_link_libraries(scontact2 mbtools2)
......@@ -46,10 +50,6 @@ target_link_libraries(scontact3fsr mbtools3)
target_compile_definitions(scontact3fsr PUBLIC -DDIMENSION=3 -DFRICTION_ENABLED=1)
add_library(scontact3f SHARED ${SRC})
target_compile_definitions(scontact3f PUBLIC -DDIMENSION=3 -DFRICTION_ENABLED=1 -DROTATION_ENABLED=1)
target_link_libraries(scontact3f mbtools3)
......
......@@ -25,6 +25,7 @@
#include "quadtree.h"
#include "vector.h"
#include <math.h>
#include <stdint.h>
#include <stdio.h>
#include <stdlib.h>
#include <float.h>
......@@ -332,10 +333,10 @@ static void _cross (const double *a, const double *b, double *c) {
static void build_basis(const double *a, double *b, double *c){
double max = fmax(fmax(fabs(a[0]), fabs(a[1])), fabs(a[2]));
if(fabs(a[0])==max){
b[0] = a[1];
b[1] = -a[0];
b[2] = 0;
if(1){
b[0] = -a[2];
b[1] = a[0];
b[2] = -a[1];
}
else if (fabs(a[1])==max){
b[0] = 0;
......@@ -444,9 +445,9 @@ static int contactParticleParticleSolve(Contact *c, ParticleProblem *p, double d
(p->velocity[c->o0*3+0]-p->velocity[c->o1*3+0])*t[0]+
(p->velocity[c->o0*3+1]-p->velocity[c->o1*3+1])*t[1]+
(p->velocity[c->o0*3+2]-p->velocity[c->o1*3+2])*t[2];
double vs = -c->cs -
(p->velocity[c->o0*3+0]-p->velocity[c->o1*3+0])*s[0]-
(p->velocity[c->o0*3+1]-p->velocity[c->o1*3+1])*s[1]-
double vs = c->cs +
(p->velocity[c->o0*3+0]-p->velocity[c->o1*3+0])*s[0]+
(p->velocity[c->o0*3+1]-p->velocity[c->o1*3+1])*s[1]+
(p->velocity[c->o0*3+2]-p->velocity[c->o1*3+2])*s[2];
#if ROTATION_ENABLED
double res0[3], res1[3];
......@@ -490,19 +491,17 @@ static int contactParticleParticleSolve(Contact *c, ParticleProblem *p, double d
#else
cs -= c->cs;
//printf("contact : %d\n",c->o0);
if(c->o0 == 827){
printf("cs : %f ct : %f dv : %f \n",cs,ct,dp);}
coordAdd(&p->velocity[c->o0 * DIMENSION], -ct*(p->a/(p->a+1))*c->a0, t);
coordAdd(&p->velocity[c->o1 * DIMENSION], ct*(p->a/(p->a+1))*c->a1, t);
coordAdd(&p->velocity[c->o0 * DIMENSION], cs*(p->a/(p->a+1))*c->a0, s);
coordAdd(&p->velocity[c->o1 * DIMENSION], -cs*(p->a/(p->a+1))*c->a1, s);
coordAdd(&p->velocity[c->o0 * DIMENSION], -cs*(p->a/(p->a+1))*c->a0, s);
coordAdd(&p->velocity[c->o1 * DIMENSION], cs*(p->a/(p->a+1))*c->a1, s);
for (int i = 0; i<3; i++){
p->omega[c->o0*3 + i] -= (1/p->particles[c->o0].r)*(1/((p->a+1)))*ct*s[i]
- (1/p->particles[c->o0].r)*(1/((p->a+1)))*cs*t[i];
p->omega[c->o1*3 + i] -= (1/p->particles[c->o1].r)*(1/((p->a+1)))*ct*s[i]
- (1/p->particles[c->o1].r)*(1/((p->a+1)))*cs*t[i];
p->omega[c->o0*3 + i] += -(1/p->particles[c->o0].r)*(1/((p->a+1)))*ct*s[i]
+ (1/p->particles[c->o0].r)*(1/((p->a+1)))*cs*t[i];
p->omega[c->o1*3 + i] += -(1/p->particles[c->o1].r)*(1/((p->a+1)))*ct*s[i]
+ (1/p->particles[c->o1].r)*(1/((p->a+1)))*cs*t[i];
}
c->cs += cs;
#endif
#else
......@@ -569,9 +568,9 @@ static int contactParticleDiskSolve(Contact *c, ParticleProblem *p, double dt, d
p->velocity[c->o1*3+0]*t[0]+
p->velocity[c->o1*3+1]*t[1]+
p->velocity[c->o1*3+2]*t[2];
double vs = -c->cs -
p->velocity[c->o1*3+0]*s[0]-
p->velocity[c->o1*3+1]*s[1]-
double vs = c->cs +
p->velocity[c->o1*3+0]*s[0]+
p->velocity[c->o1*3+1]*s[1]+
p->velocity[c->o1*3+2]*s[2];
#if ROTATION_ENABLED
......@@ -610,11 +609,12 @@ if (fabs(vt)> ((p->a+1)/p->a)*mu*dp) vt *= mu*dp/fabs(vt);
#else
cs -= c->cs;
coordAdd(&p->velocity[c->o1 * DIMENSION], -ct*(p->a/(p->a+1))*c->a1, t);
coordAdd(&p->velocity[c->o1 * DIMENSION], cs*(p->a/(p->a+1))*c->a1, s);
coordAdd(&p->velocity[c->o1 * DIMENSION], -cs*(p->a/(p->a+1))*c->a1, s);
for (int i = 0; i<3; i++){
p->omega[c->o1*3 + i] -= (1/p->particles[c->o1].r)*(1/((p->a+1)))*ct*s[i]
- (1/p->particles[c->o1].r)*(1/((p->a+1)))*cs*t[i];
p->omega[c->o1*3 + i] += -(1/p->particles[c->o1].r)*(1/((p->a+1)))*ct*s[i]
+ (1/p->particles[c->o1].r)*(1/((p->a+1)))*cs*t[i];
}
c->cs += cs;
#endif
#else
ct -= c->ct;
......@@ -672,9 +672,9 @@ static int contactParticleSegmentSolve(Contact *c, ParticleProblem *p, double dt
p->velocity[c->o1*3+0]*t[0]+
p->velocity[c->o1*3+1]*t[1]+
p->velocity[c->o1*3+2]*t[2];
double vs = -c->cs -
p->velocity[c->o1*3+0]*s[0]-
p->velocity[c->o1*3+1]*s[1]-
double vs = c->cs +
p->velocity[c->o1*3+0]*s[0]+
p->velocity[c->o1*3+1]*s[1]+
p->velocity[c->o1*3+2]*s[2];
#if ROTATION_ENABLED
......@@ -712,11 +712,12 @@ if (fabs(vt)> ((p->a+1)/p->a)*mu*dp) vt *= mu*dp/fabs(vt);
#else
cs -= c->cs;
coordAdd(&p->velocity[c->o1 * DIMENSION], -ct*(p->a/(p->a+1))*c->a1, t);
coordAdd(&p->velocity[c->o1 * DIMENSION], cs*(p->a/(p->a+1))*c->a1, s);
coordAdd(&p->velocity[c->o1 * DIMENSION], -cs*(p->a/(p->a+1))*c->a1, s);
for (int i = 0; i<3; i++){
p->omega[c->o1*3 + i] -= (1/p->particles[c->o1].r)*(1/((p->a+1)))*ct*s[i]
- (1/p->particles[c->o1].r)*(1/((p->a+1)))*cs*t[i];
p->omega[c->o1*3 + i] += -(1/p->particles[c->o1].r)*(1/((p->a+1)))*ct*s[i]
+ (1/p->particles[c->o1].r)*(1/((p->a+1)))*cs*t[i];
}
c->cs += cs;
#endif
#else
ct -= c->ct;
......@@ -762,9 +763,9 @@ static int contactParticleTriangleSolve(Contact *c, ParticleProblem *p,double dt
p->velocity[c->o1*3+1]*t[1]+
p->velocity[c->o1*3+2]*t[2];
double cs = 0;
double vs = c->cs -
p->velocity[c->o1*3+0]*s[0]-
p->velocity[c->o1*3+1]*s[1]-
double vs = c->cs +
p->velocity[c->o1*3+0]*s[0]+
p->velocity[c->o1*3+1]*s[1]+
p->velocity[c->o1*3+2]*s[2];
#if ROTATION_ENABLED
double res1[3];
......@@ -777,16 +778,12 @@ double res1[3];
(res1[0]*p->particles[c->o1].r)*s[0]+
(res1[1]*p->particles[c->o1].r)*s[1]+
(res1[2]*p->particles[c->o1].r)*s[2];
//printf("omega : %f vel : %f vs : %f \n",p->omega[3*c->o1],p->velocity[c->o1+1],vs);
#endif
const double mu = particleProblemGetMu(p,p->particleMaterial[c->o1],p->triangleMaterial[c->o0]);
if (fabs(vt)> 3.5*mu*dp) vt *= mu*dp/fabs(vt);
ct = vt;
if (fabs(vs)> 3.5*mu*dp) vs *= mu*dp/fabs(vs);
cs = vs;
//printf("vs : %f mu*dp : %f cs : %f\n",vs,mu*dp,cs);
//printf("\n");
if (dp == 0.){
cs = 0;
ct = 0;
......@@ -794,12 +791,11 @@ double res1[3];
cs -= c->cs;
ct -= c->ct;
#if ROTATION_ENABLED
//printf("cs : %f ct : %f\n",cs,ct);
coordAdd(&p->velocity[c->o1 * DIMENSION], cs*(p->a/(p->a+1))*c->a1, s);
coordAdd(&p->velocity[c->o1 * DIMENSION], -cs*(p->a/(p->a+1))*c->a1, s);
coordAdd(&p->velocity[c->o1 * DIMENSION], -ct*(p->a/(p->a+1))*c->a1, t);
for (int i = 0; i<3; i++){
p->omega[c->o1*3 + i] -= (1/p->particles[c->o1].r)*(1/((p->a+1)))*ct*s[i]
- (1/p->particles[c->o1].r)*(1/((p->a+1)))*cs*t[i];
p->omega[c->o1*3 + i] += -(1/p->particles[c->o1].r)*(1/((p->a+1)))*ct*s[i]
+ (1/p->particles[c->o1].r)*(1/((p->a+1)))*cs*t[i];
}
#else
coordAdd(&p->velocity[c->o1 * DIMENSION], -cs*c->a1, s);
......@@ -1008,7 +1004,7 @@ void ParticleToMesh(size_t nParticles, double *particles, int nElements, double
}
int listcmp(const void *i0, const void *i1) {
const u_int32_t j0 = *(u_int32_t*)i0, j1 = *(u_int32_t*)i1;
const uint32_t j0 = *(uint32_t*)i0, j1 = *(uint32_t*)i1;
return j0 - j1;
}
......@@ -1024,11 +1020,11 @@ static int intersect(const double *amin, const double *amax, const double *bmin,
static void _particleProblemGenContacts2(ParticleProblem *p, const double alert)
{
clock_t tic = clock();
u_int32_t *list = NULL;
uint32_t *list = NULL;
double bbmin[DIMENSION], bbmax[DIMENSION];
_particleProblemBoundingBox(p,bbmin,bbmax);
addAlert(alert,bbmin, bbmax);
u_int16_t N[DIMENSION];
uint16_t N[DIMENSION];
double delta = alert*6;
for (int i = 0; i < DIMENSION; ++i) {
N[i] = (bbmax[i]-bbmin[i])/delta;
......@@ -1038,7 +1034,7 @@ static void _particleProblemGenContacts2(ParticleProblem *p, const double alert)
double amin[DIMENSION], amax[DIMENSION];
particleBoundingBox(&p->particles[i], &p->position[i*DIMENSION],amin,amax);
addAlert(alert/2,amin,amax);
u_int16_t from[DIMENSION],to[DIMENSION];
uint16_t from[DIMENSION],to[DIMENSION];
for (int i = 0; i < DIMENSION; ++i){
from[i] = (amin[i]-bbmin[i])/delta;
to[i] = (amax[i]-bbmin[i])/delta;
......@@ -1046,7 +1042,7 @@ static void _particleProblemGenContacts2(ParticleProblem *p, const double alert)
#if DIMENSION == 2
for (int j = from[0]; j <= to[0]; ++j) {
for (int k = from[1]; k <= to[1]; ++k) {
u_int32_t bid = j*N[1]+k;
uint32_t bid = j*N[1]+k;
*vectorPush(&list) = bid;
*vectorPush(&list) = i;
}
......@@ -1055,7 +1051,7 @@ static void _particleProblemGenContacts2(ParticleProblem *p, const double alert)
for (int j = from[0]; j <= to[0]; ++j) {
for (int k = from[1]; k <= to[1]; ++k) {
for (int l = from[2]; l <= to[2]; ++l) {
u_int32_t bid = (j*N[1]+k)*N[2]+l;
uint32_t bid = (j*N[1]+k)*N[2]+l;
*vectorPush(&list) = bid;
*vectorPush(&list) = i;
}
......@@ -1064,7 +1060,7 @@ static void _particleProblemGenContacts2(ParticleProblem *p, const double alert)
#endif
}
clock_t tic1 = clock();
qsort(list, vectorSize(list)/2, 2*sizeof(u_int32_t), listcmp);
qsort(list, vectorSize(list)/2, 2*sizeof(uint32_t), listcmp);
clock_t tic2 = clock();
int ntot = 0;
Contact *cc=NULL;
......@@ -1173,13 +1169,13 @@ static void _particleProblemGenContacts(ParticleProblem *p, const double alert)
#else
coordAdd(&p->velocity[c->o0 * DIMENSION], -c->cs*(p->a/(p->a+1))*c->a0, s);
coordAdd(&p->velocity[c->o1 * DIMENSION], c->cs*(p->a/(p->a+1))*c->a1, s);
coordAdd(&p->velocity[c->o0 * DIMENSION], c->ct*(p->a/(p->a+1))*c->a0, t);
coordAdd(&p->velocity[c->o1 * DIMENSION], -c->ct*(p->a/(p->a+1))*c->a1, t);
coordAdd(&p->velocity[c->o0 * DIMENSION], -c->ct*(p->a/(p->a+1))*c->a0, t);
coordAdd(&p->velocity[c->o1 * DIMENSION], c->ct*(p->a/(p->a+1))*c->a1, t);
for (int i = 0; i<3; i++){
p->omega[c->o0*3 + i] -= (1/p->particles[c->o0].r)*(1/((p->a+1)))*c->ct*s[i]
- (1/p->particles[c->o0].r)*(1/((p->a+1)))*c->cs*t[i];
p->omega[c->o1*3 + i] -= (1/p->particles[c->o1].r)*(1/((p->a+1)))*c->ct*s[i]
- (1/p->particles[c->o1].r)*(1/((p->a+1)))*c->cs*t[i];
p->omega[c->o0*3 + i] += -(1/p->particles[c->o0].r)*(1/((p->a+1)))*c->ct*s[i]
+ (1/p->particles[c->o0].r)*(1/((p->a+1)))*c->cs*t[i];
p->omega[c->o1*3 + i] += -(1/p->particles[c->o1].r)*(1/((p->a+1)))*c->ct*s[i]
+ (1/p->particles[c->o1].r)*(1/((p->a+1)))*c->cs*t[i];
}
#endif
......@@ -1231,11 +1227,11 @@ static void _particleProblemGenContacts(ParticleProblem *p, const double alert)
coordAdd(&p->velocity[c->o1 * DIMENSION], c->ct*(p->a/(p->a+1))*c->a1, t);
p->omega[c->o1] -= (1/p->particles[c->o1].r)*(1/((p->a+1)))*c->ct;
#else
coordAdd(&p->velocity[c->o1 * DIMENSION], c->cs*(p->a/(p->a+1))*c->a1, s);
coordAdd(&p->velocity[c->o1 * DIMENSION], -c->cs*(p->a/(p->a+1))*c->a1, s);
coordAdd(&p->velocity[c->o1 * DIMENSION], -c->ct*(p->a/(p->a+1))*c->a1, t);
for (int i = 0; i<3; i++){
p->omega[c->o1*3 + i] -= (1/p->particles[c->o1].r)*(1/((p->a+1)))*c->ct*s[i]
- (1/p->particles[c->o1].r)*(1/((p->a+1)))*c->cs*t[i];
p->omega[c->o1*3 + i] += -(1/p->particles[c->o1].r)*(1/((p->a+1)))*c->ct*s[i]
+ (1/p->particles[c->o1].r)*(1/((p->a+1)))*c->cs*t[i];
}
#endif
#else
......@@ -1281,11 +1277,11 @@ static void _particleProblemGenContacts(ParticleProblem *p, const double alert)
coordAdd(&p->velocity[c->o1 * DIMENSION], c->ct*(p->a/(p->a+1))*c->a1, t);
p->omega[c->o1] -= (1/p->particles[c->o1].r)*(1/((p->a+1)))*c->ct;
#else
coordAdd(&p->velocity[c->o1 * DIMENSION], c->cs*(p->a/(p->a+1))*c->a1, s);
coordAdd(&p->velocity[c->o1 * DIMENSION], -c->cs*(p->a/(p->a+1))*c->a1, s);
coordAdd(&p->velocity[c->o1 * DIMENSION], -c->ct*(p->a/(p->a+1))*c->a1, t);
for (int i = 0; i<3; i++){
p->omega[c->o1*3 + i] -= (1/p->particles[c->o1].r)*(1/((p->a+1)))*c->ct*s[i]
- (1/p->particles[c->o1].r)*(1/((p->a+1)))*c->cs*t[i];
p->omega[c->o1*3 + i] += -(1/p->particles[c->o1].r)*(1/((p->a+1)))*c->ct*s[i]
+ (1/p->particles[c->o1].r)*(1/((p->a+1)))*c->cs*t[i];
}
#endif
#else
......@@ -1321,14 +1317,14 @@ static void _particleProblemGenContacts(ParticleProblem *p, const double alert)
c->cs = cold->cs;
c->ct = cold->ct;
#if ROTATION_ENABLED
coordAdd(&p->velocity[c->o1 * DIMENSION], c->cs*(p->a/(p->a+1))*c->a1, s);
coordAdd(&p->velocity[c->o1 * DIMENSION], -c->cs*(p->a/(p->a+1))*c->a1, s);
coordAdd(&p->velocity[c->o1 * DIMENSION], -c->ct*(p->a/(p->a+1))*c->a1, t);
for (int i = 0; i<3; i++){
p->omega[c->o1*3 + i] -= (1/p->particles[c->o1].r)*(1/((p->a+1)))*c->ct*s[i]
- (1/p->particles[c->o1].r)*(1/((p->a+1)))*c->cs*t[i];
p->omega[c->o1*3 + i] += -(1/p->particles[c->o1].r)*(1/((p->a+1)))*c->ct*s[i]
+ (1/p->particles[c->o1].r)*(1/((p->a+1)))*c->cs*t[i];
}
#else
coordAdd(&p->velocity[c->o1 * DIMENSION], c->cs*c->a1, s);
coordAdd(&p->velocity[c->o1 * DIMENSION], -c->cs*c->a1, s);
coordAdd(&p->velocity[c->o1 * DIMENSION], -c->ct*c->a1, t);
#endif
......@@ -1453,7 +1449,31 @@ static void _particleProblemSolveContactsQueue(ParticleProblem *p, double dt, do
Contact *c = &p->contacts[ic];
#if FRICTION_ENABLED
if(iter>4000*vectorSize(p->contacts)){ printf("Warning!\n");break;}
if(iter>2000*vectorSize(p->contacts)*(DIMENSION)){
printf("Warning : Convergence not reached\n");
for (int ik = vectorSize(p->contacts)-1; ik >= 0; --ik) {
Contact *d = &p->contacts[ik];
int conv;
switch (d->type) {
case PARTICLE_PARTICLE :
conv = contactParticleParticleSolve(d, p, dt,tol);
break;
case PARTICLE_DISK :
conv = contactParticleDiskSolve(d, p, dt, tol);
break;
case PARTICLE_SEGMENT :
conv = contactParticleSegmentSolve(d, p, dt, tol);
break;
#if DIMENSION == 3
case PARTICLE_TRIANGLE :
conv = contactParticleTriangleSolve(d, p, dt, tol);
break;
#endif
}
}
break;
}
#endif
int converged;
......
......@@ -57,6 +57,7 @@ def genInitialPosition(filename, r, rout, rhop) :
p.add_boundary_disk((0,0,0),0,"Front",material="gold")
p.add_boundary_segment((0,0,0), (0,rout,rout), "Left", material = "gold")
p.add_boundary_segment((0,0,0), (rout,rout,0), "Bottom", material = "gold")
p.add_particle((0+1.5*r,0.5,0+r), r, r**2 * np.pi * rhop,"wut");
p.write_vtk(filename,0,0)
......@@ -66,12 +67,12 @@ ii = 0
#physical parameters
g = -9.81 # gravity
rhop = 2450 # grains density
tEnd = 0.20 # final time
tEnd = 2 # final time
#numerical parameters
dt = 1e-3 # time step
#geometry parameters
rout = 1 # outer radius
rout = 10 # outer radius
r = 10e-3 # grains radius
shutil.copy("mesh.msh", outputdir +"/mesh.msh")
......@@ -85,11 +86,13 @@ p.set_friction_coefficient(5.8,"wut","gold")
#Computation loop
forces = np.zeros_like(p.velocity())
k = 0
print(sin(radians(30)))
print(sin(radians(25)))
while t < tEnd :
forces[:,2] = g*p.mass()[0]*cos(radians(25))
forces[:,0] = -g*p.mass()[0]*sin(radians(25))*(2**0.5/2)
forces[:,1] = -g*p.mass()[0]*sin(radians(25))*(2**0.5/2)
print(p.omega())
forces[:,0] = -g*p.mass()[0]*sin(radians(25))#*(2**0.5/2)
#forces[:,1] = -g*p.mass()[0]*sin(radians(25))*(2**0.5/2)
#Computation of the new velocities
vn = p.velocity() + forces * dt / p.mass()
vmax = np.max(np.sqrt(vn[:, 0]**2+ vn[:, 1]**2 + vn[:,2]**2))
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment