Skip to content

GitLab

Commit 3e28147b authored by Jonathan Lambrechts's avatar Jonathan Lambrechts
Browse files

wip

parent e388f948
Pipeline #8646 failed with stages
in 2 minutes and 33 seconds
Hide whitespace changes
Inline Side-by-side
python/scontact.py
View file @ 3e28147b
......@@ -44,11 +44,11 @@ assert(np.lib.NumpyVersion(np.__version__) >= "1.17.0")
dir_path = os.path.dirname(os.path.realpath(__file__))
lib2 = np.ctypeslib.load_library("libscontact2",dir_path)
lib2f = np.ctypeslib.load_library("libscontact2f",dir_path)
lib2fwr = np.ctypeslib.load_library("libscontact2fwr",dir_path)
lib3 = np.ctypeslib.load_library("libscontact3",dir_path)
lib3f = np.ctypeslib.load_library("libscontact3f",dir_path)
lib3fwr = np.ctypeslib.load_library("libscontact3fwr",dir_path)
#lib2f = np.ctypeslib.load_library("libscontact2f",dir_path)
#lib2fwr = np.ctypeslib.load_library("libscontact2fwr",dir_path)
#lib3 = np.ctypeslib.load_library("libscontact3",dir_path)
#lib3f = np.ctypeslib.load_library("libscontact3f",dir_path)
#lib3fwr = np.ctypeslib.load_library("libscontact3fwr",dir_path)
is_64bits = sys.maxsize > 2**32
......@@ -122,15 +122,10 @@ class ParticleProblem :
if self._p == None :
raise NameError("Cannot create particle problem.")
bndtype =[('v',np.float64,dim)]
if friction_enabled :
bndtype.append(('vt',np.float64))
if dim == 3 :
bndtype.append(('vs',np.float64))
bndtype.extend([('material',np.int32),('tag',np.int32)])
self._disktype = np.dtype(bndtype+[('x',np.float64,dim),('r',np.float64)])
self._segmenttype = np.dtype(bndtype+[('p',np.float64,(2,dim))])
self._triangletype = np.dtype(bndtype+[('p',np.float64,(3,dim))])
bndtype =[('material',np.int32),('tag',np.int32)]
self._disktype = np.dtype(bndtype+[('x',np.float64,dim),('v',np.float64,dim),('r',np.float64)])
self._segmenttype = np.dtype(bndtype+[('p',np.float64,(2,dim)),('v',np.float64,(2,dim))])
self._triangletype = np.dtype(bndtype+[('p',np.float64,(3,dim)),('v',np.float64,(3,dim))])
def __del__(self):
"""Deletes the particle solver structure."""
......@@ -254,6 +249,17 @@ class ParticleProblem :
f+= compute_fn_ft("particle_triangle",self.triangles(),tag)
return f
def set_segment_velocity(self, tag, v) :
"""Sets the velocity of a boundary to a given value.
Keyword arguments:
tag -- Name of the boundary
v -- Velocity vector to be given to the boundary
"""
self._lib.particleProblemGetTagId.restype = c_size_t
tag = self._lib.particleProblemGetTagId(self._p, tag.encode())
segments = self.segments()
segments["v"][segments["tag"]==tag] = v
def set_boundary_velocity(self, tag, v) :
"""Sets the velocity of a boundary to a given value.
Keyword arguments:
......
scontact/CMakeLists.txt
View file @ 3e28147b
......@@ -26,23 +26,23 @@ set(SRC
include_directories(. ../tools)
add_library(scontact3f SHARED ${SRC})
target_compile_definitions(scontact3f PUBLIC -DDIMENSION=3 -DFRICTION_ENABLED=1 -DROTATION_ENABLED=1)
#add_library(scontact3f SHARED ${SRC})
#target_compile_definitions(scontact3f PUBLIC -DDIMENSION=3 -DFRICTION_ENABLED=1 -DROTATION_ENABLED=1)
add_library(scontact2 SHARED ${SRC})
target_compile_definitions(scontact2 PUBLIC "-DDIMENSION=2")
add_library(scontact2f SHARED ${SRC})
target_compile_definitions(scontact2f PUBLIC -DDIMENSION=2 -DFRICTION_ENABLED=1 -DROTATION_ENABLED=1)
#add_library(scontact2f SHARED ${SRC})
#target_compile_definitions(scontact2f PUBLIC -DDIMENSION=2 -DFRICTION_ENABLED=1 -DROTATION_ENABLED=1)
add_library(scontact2fwr SHARED ${SRC})
target_compile_definitions(scontact2fwr PUBLIC -DDIMENSION=2 -DFRICTION_ENABLED=1)
#add_library(scontact2fwr SHARED ${SRC})
#target_compile_definitions(scontact2fwr PUBLIC -DDIMENSION=2 -DFRICTION_ENABLED=1)
add_library(scontact3 SHARED ${SRC})
target_compile_definitions(scontact3 PUBLIC "-DDIMENSION=3")
#add_library(scontact3 SHARED ${SRC})
#target_compile_definitions(scontact3 PUBLIC "-DDIMENSION=3")
add_library(scontact3fwr SHARED ${SRC})
target_compile_definitions(scontact3fwr PUBLIC -DDIMENSION=3 -DFRICTION_ENABLED=1)
#add_library(scontact3fwr SHARED ${SRC})
#target_compile_definitions(scontact3fwr PUBLIC -DDIMENSION=3 -DFRICTION_ENABLED=1)
......
scontact/scontact.c
View file @ 3e28147b
......@@ -155,13 +155,6 @@ static int particleInitContact(ParticleProblem *p, size_t id0, Particle *p0, dou
}
typedef struct {
double v[DIMENSION];
#if FRICTION_ENABLED
double vt;
#if DIMENSION == 3
double vs;
#endif
#endif
int material;
int tag;
}Boundary;
......@@ -169,6 +162,7 @@ typedef struct {
struct _Disk {
Boundary b;
double x[DIMENSION];
double v[DIMENSION];
double r;
};
......@@ -186,15 +180,9 @@ size_t particleProblemAddBoundaryDiskTagId(ParticleProblem *p, const double x[DI
d->b.tag = tag;
d->b.material = materialTag;
d->r = r;
#if FRICTION_ENABLED
d->b.vt = 0;
#if DIMENSION == 3
d->b.vs = 0;
#endif
#endif
for (int i = 0; i < DIMENSION; ++i) {
d->x[i] = x[i];
d->b.v[i] = 0.;
d->v[i] = 0.;
}
return vectorSize(p->disks) - 1;
}
......@@ -232,6 +220,7 @@ static int diskInitContact(size_t id, const Disk *d, size_t particle, const Part
struct _Segment{
Boundary b;
double p[2][DIMENSION];
double v[2][DIMENSION];
};
static void coordAdd(double *x, double a, const double *y) {
......@@ -249,15 +238,23 @@ static void segmentBoundingBox(const Segment *s, double *pmin, double *pmax) {
}
}
static int segmentProjectionIsInside(const void *ps, const double *x) {
static double segmentProjectionIsInside(const void *ps, const double *x, const double *n, double vnfree, double dv) {
const Segment *s = (const Segment*)ps;
double alpha = 0, beta = 0;
for (int i = 0; i < DIMENSION; ++i) {
const double d = s->p[0][i] - s->p[1][i];
const double d = s->p[1][i] - s->p[0][i];
alpha += (x[i] - s->p[0][i]) * d;
beta += (x[i] - s->p[1][i]) * d;
}
return (alpha < 0 && beta > 0);
if (alpha < 0 || beta > 0) return 0;
double vs[DIMENSION];
double xi = alpha/(alpha-beta);
for (int i = 0; i < DIMENSION; ++i) {
vs[i] = (1-xi)*s->v[0][i] + xi*s->v[1][i];
vnfree -= vs[i] * n[i];
}
if (dv > fabs(vnfree)) dv = fabs(vnfree);
return dv;
}
static int segmentInitContact(size_t id, const Segment *s, size_t particle, const Particle *p, double *x, double alert, Contact *c) {
......@@ -287,11 +284,12 @@ static int segmentInitContact(size_t id, const Segment *s, size_t particle, cons
}
contactBuildBasis(c);
c->D = nnorm - p->r;
if (c->D < 0 && segmentProjectionIsInside(s, x)) c->D = 0;
//if (c->D < 0 && segmentProjectionIsInside(s, x)) c->D = 0;
c->a0 = 0;
c->a1 = 1;
c->type = PARTICLE_SEGMENT;
return c->D >= 0 && c->D < alert;
return c->D < alert;
//return c->D >= 0 && c->D < alert;
}
......@@ -299,6 +297,7 @@ static int segmentInitContact(size_t id, const Segment *s, size_t particle, cons
struct _Triangle {
Boundary b;
double p[3][DIMENSION];
double v[3][DIMENSION];
};
void particleProblemAddBoundaryTriangleTagId(ParticleProblem *p, const double x0[3], const double x1[3], const double x2[3], int tag, int materialTag) {
......@@ -309,12 +308,10 @@ void particleProblemAddBoundaryTriangleTagId(ParticleProblem *p, const double x0
t->p[0][i] = x0[i];
t->p[1][i] = x1[i];
t->p[2][i] = x2[i];
t->b.v[i] = 0.;
t->v[0][i] = 0.;
t->v[1][i] = 0.;
t->v[2][i] = 0.;
}
#if FRICTION_ENABLED
t->b.vt = 0;
t->b.vs = 0;
#endif
}
void particleProblemAddBoundaryTriangle(ParticleProblem *p, const double x0[3], const double x1[3], const double x2[3], const char *tag, const char *material) {
......@@ -486,12 +483,6 @@ static int contactParticleParticleSolve(Contact *c, ParticleProblem *p, double d
return (fabs(dp) <= tol/dt && rfriction);
}
static double contactParticleBndSolve(Contact *c, const double *v, double vlocfree[DIMENSION], double *vn, double dt) {
for (int i = 0; i < DIMENSION; ++i)
vlocfree[i] = v[i] + c->n[i] * c->dv;
*vn = dot(vlocfree, c->n);
return fmax(0, *vn - c->D/dt);
}
#if FRICTION_ENABLED
static int contactParticleBndFrictionSolve(ParticleProblem *p, Contact *c, Boundary *bnd, double dp, double toldt) {
#if DIMENSION == 2
......@@ -562,31 +553,33 @@ static int contactParticleBndFrictionSolve(ParticleProblem *p, Contact *c, Bound
}
#endif
typedef int IsInsideFunc(const void*, const double *x);
typedef double IsInsideFunc(const void*, const double *x, const double *n, double vnfree, double dv);
static int contactParticleBoundarySolve(Contact *c, ParticleProblem *p, double dt, double tol, Boundary *bnd, IsInsideFunc *isInsideFunc) {
if(p->ForcedFlag[c->o1])
return 1;
double vlocfree[DIMENSION], vn, vloc[DIMENSION];
for (int i= 0; i < DIMENSION; ++i)
vloc[i] = p->velocity[c->o1*DIMENSION + i] - bnd->v[i];
double dp = contactParticleBndSolve(c, vloc, vlocfree, &vn, dt);
double vlocfree[DIMENSION];
double vn = 0;
for (int i= 0; i < DIMENSION; ++i) {
vlocfree[i] = p->velocity[c->o1*DIMENSION + i] + c->n[i] * c->dv;
vn += vlocfree[i] * c->n[i];
}
double dv = fmax(0, vn - c->D/dt);
double r = p->particles[c->o1].r;
if (isInsideFunc) {
double xn[DIMENSION];
for (int i = 0; i < DIMENSION; ++i)
xn[i] = p->position[c->o1*DIMENSION + i] + r*c->n[i] + vlocfree[i]*c->D / vn;
if (!isInsideFunc(bnd, xn))
dp = 0;
xn[i] = p->position[c->o1*DIMENSION + i] + r*c->n[i] + vlocfree[i]*dt;//vlocfree[i]*c->D / fmax(vn,1e-16);
dv = isInsideFunc(bnd, xn, c->n, vn, dv);
}
int rfriction = 1;
#if FRICTION_ENABLED
rfriction = contactParticleBndFrictionSolve(p,c,bnd,dp,tol/dt);
rfriction = contactParticleBndFrictionSolve(p,c,bnd,dv,tol/dt);
#endif
dp -= c->dv;
coordAdd(&p->velocity[c->o1*DIMENSION], -dp, c->n);
c->dv += dp;
return (fabs(dp) <= tol/dt && rfriction);
dv -= c->dv;
coordAdd(&p->velocity[c->o1*DIMENSION], -dv, c->n);
c->dv += dv;
return (fabs(dv) <= tol/dt && rfriction);
}
......@@ -1130,13 +1123,33 @@ int particleProblemIterate(ParticleProblem *p, double alert, double dt, double t
#endif
}
for (size_t i = 0; i < vectorSize(p->disks); ++i){
for (int j = 0; j < DIMENSION; ++j) {
p->disks[i].x[j] += p->disks[i].v[j] * dt;
}
}
for (size_t i = 0; i < vectorSize(p->segments); ++i) {
for (int j = 0; j < DIMENSION; ++j) {
p->segments[i].p[0][j] += p->segments[i].v[0][j] * dt;
p->segments[i].p[1][j] += p->segments[i].v[1][j] * dt;
}
}
#if DIMENSION == 3
for (size_t i = 0; i < vectorSize(p->triangles); ++i) {
for (int j = 0; j < DIMENSION; ++j) {
p->triangles[i].p[0][j] += p->triangles[i].v[0][j] * dt;
p->triangles[i].p[1][j] += p->triangles[i].v[1][j] * dt;
p->triangles[i].p[2][j] += p->triangles[i].v[2][j] * dt;
}
}
#endif
double *oldVelocity = malloc(sizeof(double)*vectorSize(p->velocity));
for (size_t j = 0; j < vectorSize(p->particles); ++j){
for (size_t i = 0; i < DIMENSION; ++i){
if(!p->ForcedFlag[j]){
oldVelocity[j*DIMENSION + i] = p->velocity[j * DIMENSION + i];
p->velocity[j * DIMENSION + i] += externalForces[j * DIMENSION + i] * dt / p->particles[j].m;
}
p->contactForces[j*DIMENSION+i] = p->velocity[j*DIMENSION+i];
}
......@@ -1150,37 +1163,21 @@ int particleProblemIterate(ParticleProblem *p, double alert, double dt, double t
p->contactForces[j*DIMENSION+i] = 0;
}
}
for (size_t j = 0; j<vectorSize(p->contacts); ++j){
p->contacts[j].dv /= dt;
p->contacts[j].ct /= dt;
for (size_t j = 0; j<vectorSize(p->contacts); ++j){
p->contacts[j].dv /= dt;
p->contacts[j].ct /= dt;
#if DIMENSION==3
p->contacts[j].cs /= dt;
p->contacts[j].cs /= dt;
#endif
}
}
free(oldVelocity);
return 0;
}
for (size_t i = 0; i < vectorSize(p->position); ++i){
p->position[i] += p->velocity[i] * dt;
}
for (size_t i = 0; i < vectorSize(p->disks); ++i)
for (int j = 0; j < DIMENSION; ++j)
p->disks[i].x[j] += p->disks[i].b.v[j] * dt;
for (size_t i = 0; i < vectorSize(p->segments); ++i)
for (int j = 0; j < DIMENSION; ++j) {
p->segments[i].p[0][j] += p->segments[i].b.v[j] * dt;
p->segments[i].p[1][j] += p->segments[i].b.v[j] * dt;
}
#if DIMENSION == 3
for (size_t i = 0; i < vectorSize(p->triangles); ++i) {
for (int j = 0; j < DIMENSION; ++j) {
p->triangles[i].p[0][j] += p->triangles[i].b.v[j] * dt;
p->triangles[i].p[1][j] += p->triangles[i].b.v[j] * dt;
p->triangles[i].p[2][j] += p->triangles[i].b.v[j] * dt;
}
}
#endif
for (size_t j = 0; j < vectorSize(p->particles); ++j){
for (size_t i = 0; i < DIMENSION; ++i){
p->contactForces[j*DIMENSION+i] = (p->velocity[j*DIMENSION+i]-p->contactForces[j*DIMENSION+i])*p->particles[j].m/dt;
......@@ -1206,14 +1203,9 @@ size_t particleProblemAddBoundarySegmentTagId(ParticleProblem *p, const double x
for (int i = 0; i < DIMENSION; ++i) {
s->p[0][i] = x0[i];
s->p[1][i] = x1[i];
s->b.v[i] = 0.;
s->v[0][i] = 0.;
s->v[1][i] = 0.;
}
#if FRICTION_ENABLED
s->b.vt = 0.0;
#if DIMENSION == 3
s->b.vs = 0.0;
#endif
#endif
return vectorSize(p->segments) - 1;
}
......
testcases/bndrotation/test.py 0 → 100644
View file @ 3e28147b
from migflow import scontact
import os
import numpy as np
os.makedirs("output",exist_ok=True)
p = scontact.ParticleProblem(2,False,False)
a = 1
r = 0.05
p.add_boundary_segment([0,0],[0,a],"bnd")
#p.add_boundary_segment([0,a],[-a/3,a],"bnd")
p.add_boundary_disk([-a/3,a],0,"bnd")
p.add_boundary_segment([-a,-a],[-a,a],"wall")
p.add_boundary_segment([-a,a],[a,a],"wall")
p.add_boundary_segment([a,a],[a,-a],"wall")
p.add_boundary_segment([a,-a],[-a,-a],"wall")
#p.add_boundary_segment([-a+3*r,a-2.1*r],[a-3*r,a-2.1*r],"wall")
for x in np.arange(-a+r,a-r,2*r):
for y in np.arange(-a+r,-r,2*r):
p.add_particle((x,y),r-np.random.random()*r/2,1)
i = 0
f = np.zeros_like(p.velocity())
f[:,1] -= 0.0002
p.write_vtk("output",0,0)
while True :
x = p.segments()['p'][p.segments()['tag']==0]
v = np.zeros_like(x)
v[:,:,0] = -x[:,:,1]*0.01
v[:,:,1] = x[:,:,0]*0.01
p.segments()['v'][p.segments()['tag']==0] = v
x = p.disks()['x'][p.disks()['tag']==0]
v = np.zeros_like(x)
v[:,0] = -x[:,1]*0.01
v[:,1] = x[:,0]*0.01
p.disks()['v'][p.disks()['tag']==0] = v
p.set_use_queue(0)
p.iterate(1,f,tol=1e-8,force_motion=1)
i += 1
print(i)
p.write_vtk("output",i,i)
Markdown is supported
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