Commit f39056ba authored by Jonathan Lambrechts's avatar Jonathan Lambrechts

update hydro-quebec 3d

parent 7d1cf1b1
......@@ -27,7 +27,7 @@
\usepackage{newunicodechar}
\usepackage{MnSymbol}
%\usepackage{MnSymbol}
\DeclareUnicodeCharacter{187}{{\big\rrangle}}
\DeclareUnicodeCharacter{171}{{\big\llangle}}
\DeclareUnicodeCharacter{8249}{{\big\langle}}
......
......@@ -8,11 +8,12 @@ from .legendre import *
import time
class FluidSolver :
def __init__(self, mesh, dt, rhop, rho, mu, g) :
def __init__(self, mesh, dt, rhop, rho, mu, g, imex=True) :
if isinstance(mesh, str):
mesh = Mesh(mesh)
self._jac = MeshJacobian(mesh)
self._pf = ParticleFluid(self._jac)
self._imex = imex
self._dt = dt
law = Brinkman(
mu = mu,
......@@ -21,6 +22,7 @@ class FluidSolver :
g = np.array(g),
rho = rho,
rhop = rhop,
D = self._jac.dim,
stab = dt)
d = DofManager(mesh)
if self._jac.dim == 2 :
......@@ -50,21 +52,29 @@ class FluidSolver :
def solve(self) :
tic = time.clock()
v0 = self._pf.velocity.copy()
self._law.stab = self._dt *0.5
self._sol[:] = self._solver.solve()
self._pf.velocity += self._pf.forces * self._dt*0.5 / self._pf.mass
self._law.stab = 0
self._sol[:] = self._solverEx.solve()
self._pf.velocity[:] = v0
if self._imex :
v0 = self._pf.velocity.copy()
self._law.stab = self._dt *0.5
self._sol[:] = self._solver.solve()
self._pf.velocity += self._pf.forces * self._dt*0.5 / self._pf.mass
self._law.stab = 0
self._sol[:] = self._solverEx.solve()
self._pf.velocity[:] = v0
else :
self._law.stab = self._dt
self._sol[:] = self._solver.solve()
print("cpu : ", time.clock() - tic)
return self.forces
def write_position(self, odir, oiter, time, filetype="msh"):
if filetype == "msh" :
basename = "%s/f_%05i" % (odir, oiter)
X = np.zeros((9, self._jac.x.shape[0]))
X[[0,1,3,4,6,7], :] = self._jac.x.reshape([-1, 6]).T
if self._jac.dim == 2 :
X = np.zeros(((self._jac.dim + 1) *3, self._jac.x.shape[0]))
X[[0,1,3,4,6,7], :] = self._jac.x.reshape([-1, 6]).T
else :
X = self._jac.x.reshape(-1, 12).T
self._jac.writeScalar(basename + "-xmesh.msh", "x", oiter, time, X)
else :
print("unknown export format : \"%s\"." % filetype)
......@@ -73,10 +83,10 @@ class FluidSolver :
if filetype == "msh" :
basename = "%s/f_%05i" % (odir, oiter)
self._jac.writeScalar(basename + "-porosity.msh", "porosity", oiter, time, self._jac.dofManager.getField(0, self._pf.porosity))
self._jac.writeScalar(basename + "-pressure.msh", "pressure", oiter, time, self._dof.getField(2, self._sol))
X = np.zeros((9, self._jac.x.shape[0]))
X[0::3, :] = self._dof.getField(0, self._sol)[:-1,:]
X[1::3, :] = self._dof.getField(1, self._sol)[:-1,:]
self._jac.writeScalar(basename + "-pressure.msh", "pressure", oiter, time, self._dof.getField(self._jac.dim, self._sol))
X = np.zeros(((self._jac.dim + 1) *3, self._jac.x.shape[0]))
for d in range (self._jac.dim) :
X[d::3, :] = self._dof.getField(d, self._sol)[:-1,:]
self._jac.writeScalar(basename + "-velocity.msh", "velocity", oiter, time, X)
else :
print("unknown export format : \"%s\"." % filetype)
......
from pylmgc90 import chipy
import numpy as np
import os
class LmgcInterface(object):
"""
......@@ -26,20 +27,17 @@ class LmgcInterface(object):
- writeHeader() : write header file for plot
"""
def __init__(self,space_dim, dt, theta) :
def __init__(self,space_dim, dt) :
"""
Initialize LMGC90
"""
chipy.checkDirectories()
### Set dimension in chipy for dummies ###
chipy.SetDimension(2)
chipy.Initialize()
### definition des parametres du calcul ###
chipy.utilities_logMes('INIT TIME STEPPING')
chipy.TimeEvolution_SetTimeStep(dt)
chipy.Integrator_InitTheta(theta)
### lecture du modele ###
chipy.utilities_logMes('READ BEHAVIOURS')
......@@ -80,28 +78,57 @@ class LmgcInterface(object):
self._nbDisk = chipy.DISKx_GetNbDISKx()
self._d2bMap = chipy.DISKx_GetPtrDISKx2RBDY2()
self._position = np.empty([self._nbDisk,3], 'd')
self._velocity = np.empty([self._nbDisk,3], 'd')
self._externalF= np.empty([self._nbDisk,3], 'd')
self._volume = np.empty([self._nbDisk,1], 'd')
self._mass = np.empty([self._nbDisk,1], 'd')
self._position = np.zeros([self._nbDisk,3], 'd')
self._velocity = np.zeros([self._nbDisk,3], 'd')
self._externalF= np.zeros([self._nbDisk,3], 'd')
self._volume = np.zeros([self._nbDisk,1], 'd')
self._mass = np.zeros([self._nbDisk,1], 'd')
for i in range(self._nbDisk):
self._volume[i] = np.pi * chipy.DISKx_GetContactorRadius(i+1)**2 * 2./3
self._mass[i] = chipy.RBDY2_GetBodyMass(int(self._d2bMap[i,0]))
self._tag2bnd = {}
self._tag2id = {}
for i in range(chipy.RBDY2_GetNbRBDY2()):
for t in range(chipy.RBDY2_GetNbContactor(i+1)):
tag = chipy.RBDY2_GetContactorColor(i+1,1).rstrip('_')
self._tag2id.setdefault(tag, []).append(i)
self._freq_detect = 1
self._solver_params = { 'type' :'Stored_Delassus_Loops ',
'norm' : 'Quad ',
'conv' : 1e-5,
'relax' : 1.,
'gsit1' : 200,
'gsit2' : 200
}
def setBoundaryVelocity(self, tag, v) :
print("setBoundaryVelocity not implemented")
self._tag2bnd[tag] = v
def write(self, odir, i, t) :
print("write not implemented")
chipy.WriteOutDof(1)
chipy.WriteOutVlocRloc(1)
chipy.WritePostproFiles()
chipy.WriteDisplayFiles(1, 0.01)
def iterate(self,freq_detect,freq_write,freq_display,ref_radius,**solver_args):
#def iterate(self,freq_detect,freq_write,freq_display,ref_radius,**solver_args):
def iterate(self, dt, forces):
"""
Do one step of a LMGC90 computation.
Last argument is a dictionnary holding nlgs_solver list of arguments.
"""
#
chipy.TimeEvolution_SetTimeStep(dt)
chipy.Integrator_InitTheta(1.)
self._externalF[:,:2] = forces
chipy.IncrementStep()
for tag, v in self._tag2bnd.iteritems() :
for i in self._tag2id[tag] :
for j, iv in enumerate(v):
chipy.RBDY2_SetVlocyDrivenDof(i+1,j+1,iv)
chipy.ComputeFext()
for i in range(self._nbDisk):
......@@ -109,26 +136,14 @@ class LmgcInterface(object):
chipy.ComputeBulk()
chipy.ComputeFreeVelocity()
chipy.SelectProxTactors(freq_detect)
chipy.SelectProxTactors(self._freq_detect)
chipy.RecupRloc()
chipy.ExSolver(**solver_args)
chipy.ExSolver(**self._solver_params)
chipy.StockRloc()
chipy.ComputeDof()
if freq_write > 0:
chipy.WriteOutDof(freq_write)
chipy.WriteOutVlocRloc(freq_write)
else:
chipy.WriteLastDof()
chipy.WriteLastVlocRloc()
chipy.WritePostproFiles()
if freq_display > 0:
chipy.WriteDisplayFiles(freq_display,ref_radius)
chipy.UpdateStep()
def volume(self):
......@@ -167,10 +182,80 @@ class LmgcInterface(object):
chipy.CloseDisplayFiles()
chipy.Finalize()
def scontactToLmgc90(filename):
import scontact2
from pylmgc90 import pre_lmgc
datbox_path = 'DATBOX'
if not os.path.isdir(datbox_path):
os.mkdir(datbox_path)
sc = scontact2.ParticleProblem()
sc.load(filename)
space_dim = 2
fric = 0
x = sc.position()
v = sc.velocity()
v[:] = 0.
r = sc.particles()[:, 0]
rho_particle = np.mean(sc.particles()[:, 1]/(np.pi * r**2))
avs = pre_lmgc.avatars()
mat = pre_lmgc.materials()
mod = pre_lmgc.models()
svs = pre_lmgc.see_tables()
tac = pre_lmgc.tact_behavs()
mater = pre_lmgc.material(name='STONE',type='RIGID',density=rho_particle)
model = pre_lmgc.model(name='rigid', type='MECAx', element='Rxx2D',dimension=space_dim)
mat.addMaterial(mater)
mod.addModel(model)
for i in range(r.size) :
P = pre_lmgc.rigidDisk( r=r[i], center=x[i], model=model, material=mater, color='INxxx')
P.imposeInitValue(component=[1,2],value=list(v[i,:]))
avs.addAvatar(P)
seg = sc.segments()
bndtags = set()
for i in range(seg.shape[0]) :
x0 = np.array(seg[i, 0:space_dim])
x1 = np.array(seg[i, space_dim:2*space_dim])
tag = sc.getTagName(int(sc.segmentTag()[i]))
if len(tag) > 5 :
print("maximum tag name length in lmgc is 5, '%s' is too long" % tag)
exit()
tag += "_" * (5 - len(tag))
bndtags.add(tag)
t = (x0 - x1)/np.linalg.norm(x0 - x1)
n = np.array([-t[1], t[0]]) * 1e-12
vs = np.array([x0-n, x0 + n, x1 + n, x1 - n])
av = pre_lmgc.rigidPolygon(model,mater,np.zeros([space_dim]),color=tag,generation_type="full",vertices=vs)
av.imposeDrivenDof(component=[1,2,3],dofty='vlocy')
avs.addAvatar(av)
clb_fric = pre_lmgc.tact_behav(name='iqsc0',type='IQS_CLB',fric=fric)
tac += clb_fric
svs += pre_lmgc.see_table(CorpsCandidat ="RBDY2", candidat ="DISKx", colorCandidat ='INxxx',
CorpsAntagoniste="RBDY2", antagoniste="DISKx", colorAntagoniste='INxxx',
behav=clb_fric, alert=r.min())
for tag in bndtags :
svs += pre_lmgc.see_table(CorpsCandidat ="RBDY2", candidat ="DISKx", colorCandidat ='INxxx',
CorpsAntagoniste="RBDY2", antagoniste="POLYG", colorAntagoniste=tag,
behav=clb_fric, alert=r.min())
post = pre_lmgc.postpro_commands()
solv = pre_lmgc.postpro_command(type='SOLVER INFORMATIONS', step=1)
post.addCommand(solv)
# file writting
pre_lmgc.writeBodies(avs,chemin=datbox_path)
pre_lmgc.writeDrvDof(avs,chemin=datbox_path)
pre_lmgc.writeDofIni(avs,chemin=datbox_path)
pre_lmgc.writeModels(mod,chemin=datbox_path)
pre_lmgc.writeBulkBehav(mat,chemin=datbox_path)
pre_lmgc.writeTactBehav(tac,svs,chemin=datbox_path)
pre_lmgc.writeVlocRlocIni(chemin=datbox_path)
pre_lmgc.writePostpro(commands=post, parts=avs, path=datbox_path)
......@@ -28,6 +28,7 @@ class ScontactInterface :
def setBoundaryVelocity(self, tag, v) :
p = self._p
tag = p.getTagId(tag)
p.disks()[p.diskTag() == tag, 2:4] = v
p.segments()[p.segmentTag() == tag, 4:6] = v
......
import numpy as np
import shutil
import scontact3
class ScontactInterface :
def __init__(self, inputfile) :
p = scontact3.ParticleProblem()
self._p = p
p.load(inputfile)
p.velocity()[:, :] = 0;
self._r = np.mean(p.particles()[:, 0])
self._volume = np.pi * p.particles()[:, [0]]**3 * 4./3
self._mass = p.particles()[:, [1]]
self._position = p.position()
self._velocity = p.velocity()
def volume(self):
return self._volume
def mass(self):
return self._mass
def position(self):
return self._position
def velocity(self):
return self._velocity
def setBoundaryVelocity(self, tag, v) :
p = self._p
tag = p.getTagId(tag)
p.disks()[p.diskTag() == tag, 3:6] = v
p.segments()[p.segmentTag() == tag, 6:9] = v
p.triangles()[p.triangleTag() == tag, 9:12] = v
def iterate(self, dt, forces) :
self._p.externalForces()[:] = forces
self._p.iterate(self._r, dt, 1e-5)
def write(self, odir, i, t) :
outputname = "%s/part-%05i" % (odir, i)
self._p.write(outputname + "_tmp")
shutil.move(outputname + "_tmp", outputname)
......@@ -5,6 +5,7 @@
#include <stdio.h>
#include <stdlib.h>
#include <float.h>
#include <string.h>
typedef struct _Particle Particle;
typedef struct _Contact Contact;
......@@ -19,6 +20,7 @@ struct _ParticleProblem{
Contact *contacts;
double *position, *velocity, *externalForces;
Disk *disks;
char **tagname;
int *diskTag, *segmentTag;
Segment *segments;
Contact *contactParticleDisks, *contactParticleSegments;
......@@ -100,7 +102,7 @@ static void diskBoundingBox(const Disk *d, double *pmin, double *pmax) {
}
}
size_t particleProblemAddBoundaryDisk(ParticleProblem *p, const double x[DIMENSION], double r, int tag) {
static size_t particleProblemAddBoundaryDiskTagId(ParticleProblem *p, const double x[DIMENSION], double r, size_t tag) {
Disk *d = vectorPush(&p->disks);
*vectorPush(&p->diskTag) = tag;
d->r = r;
......@@ -111,6 +113,10 @@ size_t particleProblemAddBoundaryDisk(ParticleProblem *p, const double x[DIMENSI
return vectorSize(p->disks) - 1;
}
size_t particleProblemAddBoundaryDisk(ParticleProblem *p, const double x[DIMENSION], double r, const char *tag) {
return particleProblemAddBoundaryDiskTagId(p, x, r, particleProblemGetTagId(p, tag));
}
static int diskInitContact(size_t id, const Disk *d, size_t particle, const Particle *p, double *x, double alert, Contact *c) {
double nnorm = 0;
c->dv = 0;
......@@ -187,7 +193,7 @@ struct _Triangle {
double v[DIMENSION];
};
void particleProblemAddBoundaryTriangle(ParticleProblem *p, const double x0[3], const double x1[3], const double x2[3], int tag) {
static void particleProblemAddBoundaryTriangleTagId(ParticleProblem *p, const double x0[3], const double x1[3], const double x2[3], size_t tag) {
Triangle *t = vectorPush(&p->triangles);
*vectorPush(&p->triangleTag) = tag;
for (int i = 0; i < DIMENSION; ++i) {
......@@ -198,6 +204,10 @@ void particleProblemAddBoundaryTriangle(ParticleProblem *p, const double x0[3],
}
}
void particleProblemAddBoundaryTriangle(ParticleProblem *p, const double x0[3], const double x1[3], const double x2[3], const char *tag) {
return particleProblemAddBoundaryTriangleTagId(p, x0, x1, x2, particleProblemGetTagId(p, tag));
}
static void triangleBoundingBox(const Triangle *t, double *pmin, double *pmax) {
for (int i = 0; i < 3; ++i) {
pmin[i] = fmin(fmin(t->p[0][i], t->p[1][i]), t->p[2][i]);
......@@ -688,7 +698,7 @@ void particleProblemIterate(ParticleProblem *p, double alert, double dt, double
#endif
}
size_t particleProblemAddBoundarySegment(ParticleProblem *p, const double x0[DIMENSION], const double x1[DIMENSION], int tag) {
static size_t particleProblemAddBoundarySegmentTagId(ParticleProblem *p, const double x0[DIMENSION], const double x1[DIMENSION], size_t tag) {
Segment *s = vectorPush(&p->segments);
*vectorPush(&p->segmentTag) = tag;
for (int i = 0; i < DIMENSION; ++i) {
......@@ -699,6 +709,10 @@ size_t particleProblemAddBoundarySegment(ParticleProblem *p, const double x0[DIM
return vectorSize(p->segments) - 1;
}
size_t particleProblemAddBoundarySegment(ParticleProblem *p, const double x0[DIMENSION], const double x1[DIMENSION], const char *tag) {
return particleProblemAddBoundarySegmentTagId(p, x0, x1, particleProblemGetTagId(p, tag));
}
void particleProblemAddParticle(ParticleProblem *p, const double x[DIMENSION], double r, double m) {
size_t n = vectorSize(p->particles);
Particle *particle = vectorPush(&p->particles);
......@@ -717,6 +731,11 @@ void particleProblemAddParticle(ParticleProblem *p, const double x[DIMENSION], d
void particleProblemWrite(const ParticleProblem *p, const char *filename) {
FILE *output = fopen(filename, "w");
fprintf(output, "DIMENSION %i\n", DIMENSION);
fprintf(output, "TAGS %lu", vectorSize(p->tagname));
for (size_t i = 0; i < vectorSize(p->tagname); ++i) {
fprintf(output, " %s", p->tagname[i]);
}
fprintf(output, "\n");
#if DIMENSION == 3
for (size_t i = 0; i < vectorSize(p->triangles); ++i) {
fprintf(output, "T");
......@@ -768,6 +787,7 @@ void particleProblemLoad(ParticleProblem *p, const char *filename)
FILE *input = fopen(filename, "r");
char type[128];
int tag;
char stag[256];
while (fscanf(input, "%127s", type) == 1) {
if (!strcmp(type, "DIMENSION")) {
int dim;
......@@ -777,6 +797,22 @@ void particleProblemLoad(ParticleProblem *p, const char *filename)
exit(1);
}
}
else if (!strcmp(type, "TAGS")) {
for (size_t i = 0; i < vectorSize(p->tagname); ++i)
free(p->tagname[i]);
vectorFree(p->tagname);
p->tagname = NULL;
int n;
fscanf(input, "%d", &n);
for (size_t i = 0; i < n; ++i) {
fscanf(input, "%255s", stag);
if (strlen(stag) == 255) {
printf("tag name cannot exceed 255 characters\n");
exit(1);
}
particleProblemGetTagId(p, stag);
}
}
else if (!strcmp(type, "P")) {
double point[DIMENSION], v[DIMENSION];
double r, m;
......@@ -792,14 +828,14 @@ void particleProblemLoad(ParticleProblem *p, const char *filename)
double r;
coordRead(input, point);
fscanf(input, "%le %d", &r, &tag);
particleProblemAddBoundaryDisk(p, point, r, tag);
particleProblemAddBoundaryDiskTagId(p, point, r, tag);
}
else if (!strcmp(type, "S")) {
double x0[DIMENSION], x1[DIMENSION];
coordRead(input, x0);
coordRead(input, x1);
fscanf(input, "%d", &tag);
particleProblemAddBoundarySegment(p, x0, x1, tag);
particleProblemAddBoundarySegmentTagId(p, x0, x1, tag);
}
#if DIMENSION == 3
else if (!strcmp(type, "T")) {
......@@ -808,7 +844,7 @@ void particleProblemLoad(ParticleProblem *p, const char *filename)
coordRead(input, x1);
coordRead(input, x2);
fscanf(input, "%d", &tag);
particleProblemAddBoundaryTriangle(p, x0, x1, x2, tag);
particleProblemAddBoundaryTriangleTagId(p, x0, x1, x2, tag);
}
#endif
else {
......@@ -849,9 +885,28 @@ void particleProblemDelete(ParticleProblem *p) {
vectorFree(p->triangleTag);
vectorFree(p->contactParticleTriangles);
#endif
for (size_t i = 0; i < vectorSize(p->tagname); ++i)
free(p->tagname[i]);
vectorFree(p->tagname);
free(p);
}
size_t particleProblemGetTagId(ParticleProblem *p, const char *tag) {
for (size_t i = 0; i < vectorSize(p->tagname); ++i)
if (!strcmp(tag, p->tagname[i]))
return i;
char *dup = malloc(sizeof(char) * strlen(tag));
strcpy(dup, tag);
*vectorPush(&p->tagname) = dup;
return vectorSize(p->tagname) - 1;
}
const char *particleProblemGetTagName(ParticleProblem *p, size_t id) {
if (id >= vectorSize(p->tagname))
return NULL;
return p->tagname[id];
}
ParticleProblem *particleProblemNew() {
ParticleProblem *p = (ParticleProblem*)malloc(sizeof(ParticleProblem));
p->jacobi = 0;
......@@ -859,6 +914,7 @@ ParticleProblem *particleProblemNew() {
p->particles = NULL;
p->contacts = NULL;
p->contactParticleDisks = NULL;
p->tagname = NULL;
p->contactParticleSegments = NULL;
p->disks = NULL;
p->diskTag = NULL;
......
......@@ -12,10 +12,12 @@ void particleProblemIterate(ParticleProblem *p, double alert, double dt, double
void particleProblemSolve(ParticleProblem *p, double alert, double dt, double tol, int maxit);
double particleProblemMaxDt(const ParticleProblem *p, double alert);
void particleProblemAddParticle(ParticleProblem *p, const double x[DIMENSION], double r, double m);
size_t particleProblemAddBoundaryDisk(ParticleProblem *p, const double x0[DIMENSION], double r, int tag);
size_t particleProblemAddBoundarySegment(ParticleProblem *p, const double x0[DIMENSION], const double x1[DIMENSION], int tag);
size_t particleProblemAddBoundaryDisk(ParticleProblem *p, const double x0[DIMENSION], double r, const char *tag);
size_t particleProblemAddBoundarySegment(ParticleProblem *p, const double x0[DIMENSION], const double x1[DIMENSION], const char *tag);
size_t particleProblemGetTagId(ParticleProblem *p, const char *tag);
const char* particleProblemGetTagName(ParticleProblem *p, size_t tag);
#if DIMENSION==3
void particleProblemAddBoundaryTriangle(ParticleProblem *p, const double x0[DIMENSION], const double x1[DIMENSION], const double x2[DIMENSION], int tag);
void particleProblemAddBoundaryTriangle(ParticleProblem *p, const double x0[DIMENSION], const double x1[DIMENSION], const double x2[DIMENSION], const char *tag);
#endif
double *particleProblemDisk(ParticleProblem *p);
double *particleProblemSegment(ParticleProblem *p);
......@@ -32,6 +34,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
void ParticleToMesh(size_t nParticles, double *particles, size_t nElements, double *elements, size_t *elid, double *xi);
#endif
......@@ -70,7 +70,7 @@
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));
ParticleToMesh(len, (double*)PyArray_DATA(p), PyArray_DIMS(e)[0], (double*)PyArray_DATA(e), (size_t*)PyArray_DATA(elid), (double*)PyArray_DATA(xi));
Py_DECREF(e);
Py_DECREF(p);
PyObject *r = PyTuple_Pack(2, (PyObject*) elid, (PyObject*) xi);
......@@ -101,10 +101,13 @@ struct ParticleProblem{};
particleProblemDelete($self->_this);
}
void addParticle(const double x[DIMENSION], double r, double m) { particleProblemAddParticle($self->_this, x, r, m); }
size_t addBoundaryDisk(const double x0[DIMENSION], double r, int tag) {return particleProblemAddBoundaryDisk($self->_this, x0, r, tag);}
size_t addBoundarySegment(const double x0[DIMENSION], const double x1[DIMENSION], int tag) {return particleProblemAddBoundarySegment($self->_this, x0, x1, tag);}
size_t getTagId(const char *tag) {return particleProblemGetTagId($self->_this, tag);}
const char* getTagName(size_t id){return particleProblemGetTagName($self->_this, id);}
size_t addBoundaryDisk(const double x0[DIMENSION], double r, const char *tag) {return particleProblemAddBoundaryDisk($self->_this, x0, r, tag);}
size_t addBoundarySegment(const double x0[DIMENSION], const double x1[DIMENSION], const char *tag) {return particleProblemAddBoundarySegment($self->_this, x0, x1, tag);}
#if DIMENSION==3
void addBoundaryTriangle(const double x0[DIMENSION], const double x1[DIMENSION], const double x2[DIMENSION], int tag) {particleProblemAddBoundaryTriangle($self->_this, x0, x1, x2, tag);}
void addBoundaryTriangle(const double x0[DIMENSION], const double x1[DIMENSION], const double x2[DIMENSION], const char *tag) {particleProblemAddBoundaryTriangle($self->_this, x0, x1, x2, tag);}
#endif
void write(const char *filename) { particleProblemWrite($self->_this, filename); }
void iterate(double alert, double dt, double tol = 0.005, int maxit = -1) { particleProblemIterate($self->_this, alert, dt, tol, maxit); }
......
......@@ -181,6 +181,12 @@ static void readCoord(std::istream &stream, std::vector<float> &coordVector, int
}
}
const std::string ObjectCollection::getTagName(size_t id) const {
if (id >= _tagname.size())
return "";
return _tagname[id];
}
int ObjectCollection::read(int step) {
if (step == -1)
step = _step;
......@@ -208,6 +214,17 @@ int ObjectCollection::read(int step) {
std::cout <<"Dimension not specified in file \""<<_filename<<"\"\n";
exit(0);
}
if (type == "TAGS") {
size_t n;
hf >> n;
_tagname.clear();
for (int i = 0; i < n; ++i) {
std::string buf;
hf >> buf;
_tagname.push_back(buf);
}
continue;
}
if (type == "P") {
readCoord(hf, coord, 1 * dim);
readCoord(hf, vel, 1 * dim);
......
......@@ -10,6 +10,7 @@ class ObjectCollection {
std::vector<bool> fixed;
std::vector<float> coord, segcoord, tricoord;
std::vector<float> circleGeom;
std::vector<std::string> _tagname;
int dim;
int _step;
std::string _filename;
......@@ -29,6 +30,7 @@ class ObjectCollection {
int read(int step = -1);
void display();
const std::string getTagName(size_t id) const;
void shift(double dx, double dy);
void zoom_at(double factor, double cx, double cy);
void rotate(double r0, double r1);
......
......@@ -171,14 +171,14 @@ static void update_menus(Fl_Menu_ *menubar, ScontactPlotWindow *wingl) {
int idx = menubar->find_index("&Display/&Visibility");
menubar->clear_submenu(idx); for (std::map<int, int>::iterator it = collection.visible_flag().begin(); it!= collection.visible_flag().end(); ++it) {
std::ostringstream oss;
oss << "&" << it->first;
oss << "&" << it->first << " " << wingl->collection().getTagName(it->first);
menubar->insert(++idx, oss.str().c_str(), 0, visibility_cb, wingl, FL_MENU_TOGGLE | ((!(it->second &2)) ? FL_MENU_VALUE : 0));
}
idx = menubar->find_index("&Display/Visibility &Boundary");
menubar->clear_submenu(idx);
for (std::map<int, int>::iterator it = collection.visible_flag().begin(); it!= collection.visible_flag().end(); ++it) {
std::ostringstream oss;
oss << "&" << it->first;
oss << "&" << it->first << " " << wingl->collection().getTagName(it->first);
menubar->insert(++idx, oss.str().c_str(), 0, visibility_bnd_cb, wingl, FL_MENU_TOGGLE | ((!it->second &1) ? FL_MENU_VALUE : 0));
}
Fl_Menu_Item *item = (Fl_Menu_Item*) menubar->find_item("&Display/&Clip");
......
import scontact2
import numpy as np
import mesh
from pyfefp import mesh
import random
import os, time, sys
import shutil, random
......@@ -11,27 +11,27 @@ if not os.path.isdir(outputdir) :