Commit cf7f249c authored by Jonathan Lambrechts's avatar Jonathan Lambrechts

rewrite fe code in python/numpy + no stabilisation + correct mass

lumping
parent 6e9bdfe5
DOCS= particle-fluid-interaction ale equations stability
DOCS= particle-fluid-interaction ale equations stability stability2
all : ${DOCS:=.pdf}
%.pdf : %.tex
rubber --pdf $<
......
......@@ -340,6 +340,102 @@ static void _particleProblemBoundingBox(ParticleProblem *p, double *bbmin, doubl
#endif
}
static void _bboxtol(double *bbmin, double *bbmax) {
double lmax = 0;
for (int i = 0; i < DIMENSION; ++i) {
lmax = fmax(lmax, bbmax[i] - bbmin[i]);
}
double tol = 1e-8 * lmax;
for (int i = 0; i < DIMENSION; ++i) {
bbmax[i] += tol;
bbmin[i] -= tol;
}
}
void ParticleToMesh(size_t nParticles, double *particles, size_t nElements, double *elements, size_t *elid, double *Xi)
{
double bbmin[DIMENSION], bbmax[DIMENSION];
int N = DIMENSION + 1;
for (int i = 0; i < DIMENSION; ++i) {
bbmin[i] = elements[i];
bbmax[i] = elements[i];
}
for (size_t i = 0; i < N * nElements; ++i) {
for (int d = 0; d < DIMENSION; ++d) {
bbmin[d] = fmin(bbmin[d], elements[DIMENSION * i + d]);
bbmax[d] = fmax(bbmax[d], elements[DIMENSION * i + d]);
}
}
_bboxtol(bbmin, bbmax);
Cell *tree = cellNew(bbmin, bbmax);
double amin[DIMENSION], amax[DIMENSION];
for (size_t i = 0; i < nElements; ++i) {
double *el = elements + (DIMENSION * N) * i;
for (int d = 0; d < DIMENSION; ++d) {
amin[d] = el[d];
amax[d] = el[d];
}
for (int v = 1; v < N; ++v) {//only simplices
for (int d = 0; d < DIMENSION; ++d) {
amin[d] = fmin(amin[d], el[v * DIMENSION + d]);
amax[d] = fmax(amax[d], el[v * DIMENSION + d]);
}
}
_bboxtol(amin, amax);
cellAdd(tree, amin, amax, i, NULL);
}
int *found = NULL;
vectorPushN(&found, 100);
vectorClear(found);
for (size_t i = 0; i < nParticles; ++i) {
double *x = &particles[i * DIMENSION];
elid[i] = -1;
cellSearch(tree, x, x, &found);
for (size_t j = 0; j < vectorSize(found); ++j) {
double *el = &elements[found[j] * N * DIMENSION];
double *X[3] = {el, el + DIMENSION, el + DIMENSION * 2};
double dx = x[0] - X[0][0], dy = x[1] - X[0][1];
double DX[2] = {X[1][0] - X[0][0], X[2][0] - X[0][0]};
double DY[2] = {X[1][1] - X[0][1], X[2][1] - X[0][1]};
double det = DX[1] * DY[0] - DX[0] * DY[1];
double xi = (DX[1] * dy - DY[1] * dx) / det;
double eta = -(DX[0] * dy - DY[0] * dx) / det;
if (xi > -1e-8 && eta > -1e-8 && 1 - xi - eta > -1e-8) {
Xi[i * DIMENSION + 0] = xi;
Xi[i * DIMENSION + 1] = eta;
elid[i] = found[j];
break;
}
}
if (elid[i] == -1) {
printf(" PARTICLE %i OUTSIDE DOMAIN N2 search!!!\n", i);
double toll = -1000;
for (size_t j = 0; j < nElements; ++j) {
double *el = &elements[j * N * DIMENSION];
double *X[3] = {el, el + DIMENSION, el + DIMENSION * 2};
double dx = x[0] - X[0][0], dy = x[1] - X[0][1];
double DX[2] = {X[1][0] - X[0][0], X[2][0] - X[0][0]};
double DY[2] = {X[1][1] - X[0][1], X[2][1] - X[0][1]};
double det = DX[1] * DY[0] - DX[0] * DY[1];
double xi = (DX[1] * dy - DY[1] * dx) / det;
double eta = -(DX[0] * dy - DY[0] * dx) / det;
toll = fmax(toll, fmin(-xi, fmin(-eta, xi + eta -1)));
if (xi > -1e-8 && eta > -1e-8 && 1 - xi - eta > -1e-8) {
Xi[i * DIMENSION + 0] = xi;
Xi[i * DIMENSION + 1] = eta;
elid[i] = j;
break;
}
}
if (elid[i] == -1) {
printf(" PARTICLE %i OUTSIDE DOMAIN!!! %g %g (%g)\n", i, x[0], x[1], toll);
exit(1);
}
}
}
cellFree(tree);
}
static void _particleProblemGenContacts(ParticleProblem *p, const double alert)
{
size_t iold = 0;
......@@ -730,7 +826,3 @@ unsigned long int particleProblemNParticle(ParticleProblem *p) {
double *particleProblemRadius(ParticleProblem *p) {
return (double*)p->particles;
}
......@@ -31,5 +31,6 @@ int *particleProblemSegmentTag(ParticleProblem *p);
#if DIMENSION == 3
int *particleProblemTriangleTag(ParticleProblem *p);
double *particleProblemTriangle(ParticleProblem *p);
void ParticleToMesh(size_t nParticles, double *particles, size_t nElements, double *elements, size_t *elid, double *xi);
#endif
#endif
......@@ -51,9 +51,30 @@
}
return 1;
}
struct ParticleProblem{
ParticleProblem *_this;
};
struct ParticleProblem{
ParticleProblem *_this;
};
PyObject *particleToMesh(PyObject *position, PyObject *elements) {
PyArrayObject *p = (PyArrayObject*)PyArray_FROMANY(position, NPY_DOUBLE, 2, 2, NPY_ARRAY_CARRAY_RO);
if (!p || PyArray_DIMS(p)[1] != DIMENSION) {
PyErr_Format(PyExc_TypeError, "cannot convert argument 1 of particleToMesh to an array of coordinates of dimension %i", DIMENSION);
return NULL;
}
PyArrayObject *e = (PyArrayObject*)PyArray_FROMANY(elements, NPY_DOUBLE, 3, 3, NPY_ARRAY_CARRAY_RO);
if (!e || PyArray_DIMS(e)[1] != (DIMENSION + 1) || PyArray_DIMS(e)[2] != DIMENSION) {
Py_DECREF(p);
PyErr_Format(PyExc_TypeError, "cannot convert argument 3 of particleToMesh to an array of elements (double[][%i][%i])", DIMENSION +1 , DIMENSION);
return NULL;
}
npy_intp len = PyArray_DIMS(p)[0];
PyArrayObject *elid = (PyArrayObject*)PyArray_New(&PyArray_Type, 1, &len, NPY_UINTP, NULL, NULL, 0, 0, NULL);
npy_intp dd[2] = {len, DIMENSION};
PyArrayObject *xi = (PyArrayObject*)PyArray_New(&PyArray_Type, 2, dd, NPY_DOUBLE, NULL, NULL, 0, 0, NULL);
ParticleToMesh(len, (double*)PyArray_DATA(p), PyArray_DIMS(e)[0], (size_t*)PyArray_DATA(e), (size_t*)PyArray_DATA(elid), (double*)PyArray_DATA(xi));
Py_DECREF(e);
Py_DECREF(p);
return PyTuple_Pack(2, (PyObject*) elid, (PyObject*) xi);
}
%}
%typemap(in) const double[DIMENSION] (double tmp[DIMENSION]){
......@@ -113,3 +134,5 @@ struct ParticleProblem{};
void useJacobi(int jacobi, double relax = 0.5) {particleProblemUseJacobi($self->_this, jacobi, relax);}
void writeVtu(const char *filename) {particleProblemWriteVtu($self->_this, filename);}
};
PyObject *particleToMesh(PyObject *position, PyObject *elements);
......@@ -109,6 +109,7 @@ class ObjectCollection {
}
int read(const std::string basename, int step)
{
printf("%i\n", step);
std::string filename;
std::ostringstream oss;
oss << basename << "/part-" <<std::setw(5)<<std::setfill('0')<< step;
......@@ -151,10 +152,14 @@ class ObjectCollection {
else if (type == "S"){
readCoord(hf, segcoord, 2);
hf >> tag;
//if (tag != 4)
// segcoord.resize(segcoord.size() - 2 * dim);
}
else if (type == "T"){
readCoord(hf, tricoord, 3);
hf >> tag;
//if (tag != 4)
// tricoord.resize(tricoord.size() - 3 * dim);
}
else {
std::cout <<"Unknown type \""<< type << "\"in file \""<<filename<<"\"\n";
......@@ -219,7 +224,7 @@ class ObjectCollection {
}
glEnableClientState(GL_VERTEX_ARRAY);
glColor4f(0., 1., 1., 1);
glPolygonMode( GL_FRONT_AND_BACK, GL_LINE );
glPolygonMode( GL_FRONT_AND_BACK, GL_FILL ); // GL_LINE
glVertexPointer(dim, GL_FLOAT, 0, &tricoord[0]);
glDrawArrays(GL_TRIANGLES, 0, tricoord.size() / dim);
glPolygonMode( GL_FRONT_AND_BACK, GL_FILL );
......@@ -238,6 +243,51 @@ std::string basename;
FILE* ffmpeg = NULL;
ObjectCollection collection;
/* glut main loop */
#include <GL/glut.h>
int paused = 0;
int istep = -1, lstep = 0;
int step = 1;
double alpha = 0, beta = 0, scale = 1, shift[2] = {0, 0};
void idle_cb() {
if (! paused && lstep == istep) {
istep += step;
}
if (collection.read(basename, istep)) {
lstep = istep;
glutPostRedisplay();
}
}
void save_opt(const char *filename) {
FILE *f = fopen(filename, "w");
fprintf(f, "alpha %.16g\nbeta %1.6g\nscale %.16g\nshift %.16g %.16g\n", alpha, beta, scale, shift[0], shift[1]);
fclose(f);
}
void load_opt(const char *filename) {
std::ifstream f(filename);
if (!f) {
std::cerr << "warning : cannot open file \""<< filename << "\"\n";
return;
}
std::string word;
while (f>>word) {
if (word == "alpha") f >> alpha;
else if (word == "beta") f >> beta;
else if (word == "scale") f >> scale;
else if (word == "shift") f >> shift[0] >> shift[1];
else {
std::cerr << "warning : unknown option \""<< word << "\"\n";
f.ignore(30000, '\n');
}
}
}
/* command line options */
#include <argp.h>
......@@ -247,9 +297,10 @@ static char doc[] = "scontactplot -- plot the output of scontact";
static struct argp_option options[] = {
{"movie", 'm', "FILE", 0, "generate a mp4 file using ffmpeg"},
{"step", 's', "INT", 0, "step between part file"},
{"load", 'l', "FILE", 0, "load option file"},
{"load", 'i', "INT", 0, "initial step"},
{0}
};
int step = 1;
static char args_doc[] = "DIRECTORY";
static error_t parse_opt(int key, char *arg, struct argp_state *state)
{
......@@ -263,8 +314,18 @@ static error_t parse_opt(int key, char *arg, struct argp_state *state)
case 's' :{
if (sscanf(arg, "%d", &step) != 1)
return ARGP_ERR_UNKNOWN;
lstep = istep = step;
}
break;
case 'l' :{
load_opt(arg);
}
break;
case 'i' :{
if (sscanf(arg, "%d", &istep) != 1)
return ARGP_ERR_UNKNOWN;
lstep = istep;
}
case ARGP_KEY_ARG :
basename = arg;
break;
......@@ -277,26 +338,8 @@ static error_t parse_opt(int key, char *arg, struct argp_state *state)
}
return 0;
}
static struct argp argp = {options, parse_opt, args_doc, doc};
/* glut main loop */
#include <GL/glut.h>
int paused = 0;
int istep = -1, lstep = 0;
double alpha = 0, beta = 0, scale = 1, shift[2] = {0, 0};
void idle_cb() {
if (! paused && lstep == istep) {
istep += step;
}
if (collection.read(basename, istep)) {
lstep = istep;
glutPostRedisplay();
}
}
void key_cb(unsigned char key, int x, int y)
{
if (glutGetModifiers() == 0) {
......@@ -308,6 +351,9 @@ void key_cb(unsigned char key, int x, int y)
istep = 0;
glutPostRedisplay();
}
if (key == 's') {
save_opt("scontact.opt");
}
if (key == 27) {
istep = 0;
exit(0);
......@@ -368,8 +414,9 @@ void display_cb() {
}
int main(int argc, char **argv) {
argp_parse(&argp, argc, argv, 0, 0, 0);
lstep = istep = step;
scale = -1;
argp_parse(&argp, argc, argv, 0, 0, 0);
glutInit(&argc, argv);
glutInitDisplayMode(GLUT_RGBA | GLUT_DOUBLE |GLUT_DEPTH|GLUT_MULTISAMPLE);
glutInitWindowSize(1000, 1000);
......@@ -377,7 +424,8 @@ int main(int argc, char **argv) {
collection.read(basename, istep);
float bbox[6];
collection.compute_bbox(bbox);
scale = std::max(bbox[1] - bbox[0], std::max(bbox[3] - bbox[2], bbox[5] - bbox[4]))/2;
if (scale == -1)
scale = std::max(bbox[1] - bbox[0], std::max(bbox[3] - bbox[2], bbox[5] - bbox[4]))/2;
glutDisplayFunc(display_cb);
glutPostRedisplay();
glutKeyboardFunc(key_cb);
......
import numpy as np
import quadrature
D = 2
_u_ = 0
_v_ = 1
_w_ = 2
_p_ = D
_U_ = range(D)
def derivate(f, x, nd = 1, eps = 1e-8) :
f0 = f(x)
xd = x.reshape(x.shape[:x.ndim-nd]+(-1,))
d = np.ndarray((f0.shape +(xd.shape[-1],)))
for i in range(xd.shape[-1]) :
xd[..., i] += eps
d[..., i] = (f(xd.reshape(x.shape)) - f0) / eps
xd[..., i] -= eps
return d.reshape(f0.shape + x.shape[x.ndim - nd:])
def dot(a, b) :
return a[..., 0] * b[..., 0] + a[..., 1] * b[..., 1] + (0 if D == 2 else a[..., 2] * b[..., 2])
class brinkman :
def psiTerm0(self, sol, solgrad):
val = np.zeros_like(sol)
val[..., _U_] = - solgrad[..., _p_, :] * self.ca + self.mu * (\
+ dot(solgrad[..., _U_, :], self.dca) \
- sol[..., _U_] / self.ca **2 * dot(self.dca, self.dca))
val[..., [_v_]] += -self.ca
val[..., _p_] = - solgrad[..., _u_, 0] - solgrad[..., _v_, 1] - (0 if D == 2 else solgrad[..., _w_, 2])
return val
def psiTerm1(self, sol, solgrad):
val = np.ndarray(sol.shape + (D,))
val[..., _U_, :] = -solgrad[..., _U_, :] * self.mu
val[..., _p_, :] = self.vs
if self.ns :
val[..., _u_, :] += - sol[..., [_u_]] * sol[..., _U_] / self.ca
val[..., _v_, :] += - sol[..., [_v_]] * sol[..., _U_] / self.ca
if D == 3:
val[..., _w_, :] += - sol[..., [_w_]] * sol[..., _U_] / self.ca
return val
def particleTerm(self, sol, solgrad) :
val = np.zeros_like(sol)
if D == 2 :
d = 2 * (self.pf.volume[0, 0] / np.pi)**0.5
else :
d = 2 * (3./(4 * np.pi) * self.pf.volume[0, 0])**(1./3)
U = self.vs - sol[..., _U_]/self.ca
normU = np.sqrt(dot(U, U))[...,np.newaxis]
Re_O_u = d * self.ca/self.mu * self.rhof
Cd_u = (0.63 * normU + 4.8/(Re_O_u**0.5))**2
f = self.ca**(-1.8)
F = solgrad[..., _p_, :]
G = f * Cd_u * self.rhof * U / 2
val[..., _U_] = G
val[..., _p_] = 0
return val
def __init__(self, mu, pf, ns = False, rhof = 1) :
self.ns = ns
self.rhof = 1
self.mu = mu
self.pf = pf
self.strongBC = {}
def particleInteractionTerm(self, meshJ, solution, res, jac) :
dofm = solution.dof_manager()
self.ca = self.pf.porosity.at_points(self.pf.eid, self.pf.uvw)
self.vs = self.pf.velocity
sol = solution.at_points(self.pf.eid, self.pf.uvw)
gradsol = solution.gradient_at_points(self.pf.eid, self.pf.uvw, meshJ)
pTerm = self.particleTerm(sol, gradsol) * self.pf.volume
self.pf.forces[:, :] = -pTerm[:, _U_] - gradsol[..., _p_, :] * self.pf.volume - self.pf.volume * np.array([0, 1])
for i in range(dofm.n_fields()) :
np.add.at(res, dofm._fields[i][0][:, self.pf.eid], (dofm.get_field_basis(i).f(self.pf.uvw)*pTerm[:, [i]]).T)
d00 = derivate(lambda s : self.particleTerm(s, gradsol) * self.pf.volume, sol, 1)
d01 = derivate(lambda s : self.particleTerm(sol, s) * self.pf.volume, gradsol, 2)
for f in range(dofm.n_fields()) :
ff = dofm.get_field_basis(f).f(self.pf.uvw)
for g in range(dofm.n_fields()) :
pp = d00[:, f, g, None] * dofm.get_field_basis(g).f(self.pf.uvw)
if (np.max(np.abs(d01[:,f,g,:])) > 1e-8) :
d01xi = meshJ.multDXiDXLeft(d01[:, f, g, :], 1, self.pf.eid)
pp += (d01xi[:, None, :] * dofm.get_field_basis(g).df(self.pf.uvw)).sum(2)
jac.add_to_field(f, g, pp[:, None, :] * ff[:, :, None], self.pf.eid)
def fluidTerm(self, meshJ, solution, res, jac) :
dofm = solution.dof_manager()
qp, qw = quadrature.get_triangle(3)
self.vs = self.pf.momentum.at_qp(qp)
self.dvs = self.pf.momentum.gradient_at_qp(qp, meshJ)
self.ca = self.pf.porosity.at_qp(qp)
self.dca = self.pf.porosity.gradient_at_qp(qp, meshJ)
sol = solution.at_qp(qp)
gradsol = solution.gradient_at_qp(qp, meshJ)
x0 = meshJ.multJ(self.psiTerm0(sol, gradsol))
x1 = meshJ.multJ(meshJ.multDXiDXLeft(self.psiTerm1(sol, gradsol), 3))
for i in range(dofm.n_fields()) :
basis = dofm.get_field_basis(i)
psiW = basis.f(qp) * qw[:, None]
r0 = np.tensordot(psiW, x0[:, :, i], [0, 1])
dPsiW = basis.df(qp) * qw[:, None, None]
r1 = np.tensordot(dPsiW, x1[:, :, i, :], [[0, 2], [1, 2]])
np.add.at(res, dofm._fields[i][0], r0 + r1)
d00 = derivate(lambda s : self.psiTerm0(s, gradsol), sol, 1)
d01 = derivate(lambda s : self.psiTerm0(sol, s), gradsol, 2)
d10 = derivate(lambda s : self.psiTerm1(s, gradsol) , sol, 1)
d11 = derivate(lambda s : self.psiTerm1(sol, s), gradsol, 2)
for f in range(dofm.n_fields()):
basisf = dofm.get_field_basis(f)
for g in range(dofm.n_fields()):
basisg = dofm.get_field_basis(g)
if (np.max(np.abs(d00[...,f,g])) > 1e-8) :
psiPsiW = np.einsum("qn, qm, q-> qnm", basisf.f(qp), basisg.f(qp), qw)
jac.add_to_field(f, g, np.tensordot(meshJ.multJ(d00[:,:,f,g]), psiPsiW, 1))
if (np.max(np.abs(d01[:,:,f,g,:])) > 1e-8) :
dPsiPsiW = np.einsum("qmd, qn, q -> qdnm", basisg.df(qp), basisf.f(qp), qw)
jac.add_to_field(f, g, np.tensordot(meshJ.multJ(meshJ.multDXiDXLeft(d01[:,:,f,g,:], 2)), dPsiPsiW, 2))
if (np.max(np.abs(d10[:,:,f,:,g])) > 1e-8) :
psiDPsiW = np.einsum("qmd, qn, q -> qdmn", basisf.df(qp), basisg.f(qp), qw)
jac.add_to_field(f, g, np.tensordot(meshJ.multJ(meshJ.multDXiDXLeft(d10[:,:,f,:,g], 2)), psiDPsiW, 2))
if (np.max(np.abs(d11[:,:,f,:,g,:])) > 1e-8) :
dPsiDPsiW = np.einsum("qmx, qny, q -> qxymn", basisf.df(qp), basisg.df(qp), qw)
jac.add_to_field(f, g, np.tensordot(meshJ.multDXiDXLeft(meshJ.multDXiDXLeft(meshJ.multJ(d11[:,:,f,:,g,:]), 3), 2), dPsiDPsiW, 3))
def volumeTerm(self, meshJ, solution) :
dofm = solution.dof_manager()
jac = dofm.new_matrix()
res = dofm.new_vector()
self.particleInteractionTerm(meshJ, solution, res, jac)
self.fluidTerm(meshJ, solution, res, jac)
return res, jac.to_csr_matrix()
import numpy as np
import time
import scipy.sparse
import scipy.sparse.linalg
class steady :
def __init__(self, law, dof, jac) :
self.law = law
self.meshJacobian = jac
self.dof = dof
self.preconditioner = None
self.tspsolve = 0
self.tvt = 0
self.tas = 0
def solve(self) :
x = self.dof.new_vector()
x[self.dof._fixedidx] = self.dof._fixedval
firstR = None
newtonIter = 0
ttot = time.clock()
while True :
tic = time.clock()
res, jac = self.law.volumeTerm(self.meshJacobian, x)
self.tvt += time.clock() - tic
tic = time.clock()
for row in self.dof._fixedidx :
a = slice(jac.indptr[row], jac.indptr[row +1])
jac.data[a] = np.where(jac.indices[a] == row, 1., 0.)
jac.eliminate_zeros()
self.tas += time.clock() - tic
res[self.dof._fixedidx] = 0
resn = np.linalg.norm(res)
if firstR :
print("newton (%i) % 6.2e %.1f" % (newtonIter, resn, -np.log10(resn/firstR)))
if (resn/firstR < 1e-6) :
break
else :
firstR = resn
print("newton (%i) % 6.2e" % (newtonIter, resn))
newtonIter += 1
if newtonIter > 10:
break
tic = time.clock()
matrix = jac
if not self.preconditioner :
P = scipy.sparse.linalg.spilu(matrix)
self.preconditioner = scipy.sparse.linalg.LinearOperator(matrix.shape, lambda x: P.solve(x))
#I think there is a bug in scipy, without the callback argument, maxiter as no effect
dx, info = scipy.sparse.linalg.gmres(matrix, res, tol=1e-10, M=self.preconditioner, maxiter=15, callback=lambda x:None)
if info != 0 :
print("*** PRECOND ***")
P = scipy.sparse.linalg.spilu(matrix)
self.preconditioner = scipy.sparse.linalg.LinearOperator(matrix.shape, lambda x: P.solve(x))
dx,info = scipy.sparse.linalg.gmres(matrix, res, tol=1e-10, M=self.preconditioner)
x -= dx
self.tspsolve += time.clock() - tic
print("spsolve: %.3g term: %.3g assemble: %.3g" % (self.tspsolve, self.tvt, self.tas))
return x
#!/usr/bin/env python
import scontact2
import numpy as np
import mesh
import random
import os, time, sys
import shutil, random
outputdir ="output"
if not os.path.isdir(outputdir) :
os.makedirs(outputdir)
addv = set()
p = scontact2.ParticleProblem()
def loadMeshBoundary(p, filename, tags, shift= [0, 0]) :
m = mesh.mesh(filename)
pnum = [m.getPhysicalNumber(1, i) for i in tags]
for ent in m.entities :
if ent.dimension != 1 or not [i for i in pnum if i in ent.physicals]:
continue
tag = 0 if not ent.physicals else ent.physicals[0]
for i, el in enumerate(ent.elements) :
for v in el.vertices :
if v[3] in addv :
continue
addv.add(v[3])
p.addBoundaryDisk((v[0] + shift[0], v[1] + shift[1]), 0.000, tag)
v0 = el.vertices[1]
v1 = el.vertices[0]
p.addBoundarySegment((v0[0] + shift[0], v0[1] + shift[1]), (v1[0] + shift[0], v1[1] + shift[1]), tag)
def genInitialPosition(p, N, rmin, rmax, L, H, rhop) :
x = np.arange(L - rmax, -L + rmax, 2 * -rmax)
y = -np.arange(-H+rmin, H, 2* rmax);
x, y = np.meshgrid(x, y)
print(x.shape)
R = rmin + (rmax - rmin ) * np.random.random([N])
x = x.flat[:N]
y = y.flat[:N]
for i in range(N) :
p.addParticle((x[i], y[i]), R[i], R[i]**2 * np.pi * rhop)
return p
r = 0.013*2
dr = 0#0.003*2
N = 2000
genInitialPosition(p, N = N, rmin = r -dr, rmax = r + dr , L = 2, H = 1.5, rhop=2)
p.useJacobi(False)
loadMeshBoundary(p, "mesh.msh", ["Box"])
g = 1
tEnd = 10
t = 0
iter = 0
tic = time.clock()
outputname = "%s/part-%05i" % (outputdir, 0)
p.write(outputname)
ioutput = 0
p.velocity()[:, 1] = 0;
while t < tEnd :
p.externalForces()[:, :] = - p.velocity()
p.externalForces()[:, 1] += -g
p.externalForces()[:, :] *= p.particles()[:, [1]]
dt = min(0.001, p.maxdt(r +dr)/2)
p.iterate(r + dr, dt, 1e-5)
t += dt
iter += 1
if iter %50 == 0 :
outputname = "%s/part-%05i" % (outputdir, ioutput)
tmpoutputname = "%s/part-%05i_tmp" % (outputdir, ioutput)
ioutput += 1
p.write(tmpoutputname)
shutil.move(tmpoutputname, outputname)
print("cpu time : %g" % (time.clock() - tic))
loadMeshBoundary(p, "mesh.msh", ["Pile"], [0, 0.5])
p.write("deposited")
import mesh
import numpy as np
import scipy.sparse
class dofManager :