Commit 012e9ee2 authored by Jonathan Lambrechts's avatar Jonathan Lambrechts

add hydro-quebec3d

parent f7f20d4a
......@@ -184,6 +184,7 @@ static int segmentInitContact(size_t id, const Segment *s, size_t particle, cons
#if DIMENSION == 3
struct _Triangle {
double p[3][DIMENSION];
double v[DIMENSION];
};
void particleProblemAddBoundaryTriangle(ParticleProblem *p, const double x0[3], const double x1[3], const double x2[3], int tag) {
......@@ -193,6 +194,7 @@ void particleProblemAddBoundaryTriangle(ParticleProblem *p, const double x0[3],
t->p[0][i] = x0[i];
t->p[1][i] = x1[i];
t->p[2][i] = x2[i];
t->v[i] = 0.;
}
}
......@@ -253,6 +255,8 @@ static int triangleInitContact(size_t id, const Triangle *t, size_t particle, co
c->D -= p->r;
c->a0 = 0.;
c->a1 = 1.;
if (c->D < 0)
c->D = 0;
return (c->D > 0 || triangleProjectionIsInside(t, x)) && c->D < alert;
}
......@@ -432,9 +436,9 @@ static void _particleProblemGenContacts(ParticleProblem *p, const double alert)
cellFree(tree);
}
static void _particleProblemSolveContacts(ParticleProblem *p, double dt, double tol) {
static void _particleProblemSolveContacts(ParticleProblem *p, double dt, double tol, int maxit) {
int iter = 0;
for (double residual = DBL_MAX; residual > tol/dt;) {
for (double residual = DBL_MAX; residual > tol/dt && (maxit < 0 || iter < maxit);) {
residual = 0;
iter++;
double *velocityNew = p->velocity;
......@@ -454,10 +458,12 @@ static void _particleProblemSolveContacts(ParticleProblem *p, double dt, double
#if DIMENSION == 3
for (size_t ii = 0; ii <vectorSize(p->contactParticleTriangles); ++ii) {
Contact *c = &p->contactParticleTriangles[ii];
double vlocfree[DIMENSION], vn, xn[DIMENSION];
double pp = contactParticleBndSolve(c, &p->velocity[c->o1 * DIMENSION], vlocfree, &vn, dt);
double vlocfree[DIMENSION], vn, xn[DIMENSION], vloc[DIMENSION];
for (int i= 0; i < DIMENSION; ++i)
vloc[i] = p->velocity[c->o1 * DIMENSION + i] - p->triangles[c->o0].v[i];
double pp = contactParticleBndSolve(c, vloc, vlocfree, &vn, dt);
for (int i = 0; i < DIMENSION; ++i)
xn[i] = p->position[c->o1 * DIMENSION + i] - c->n[i] * p->particles[c->o1].r + vlocfree[i] * c->D/vn;
xn[i] = p->position[c->o1 * DIMENSION + i] + p->particles[c->o1].r * c->n[i] + vlocfree[i] * c->D/vn;
if (!triangleProjectionIsInside(&p->triangles[c->o0], xn))
pp = 0;
double dp = pp - c->dv;
......@@ -474,13 +480,10 @@ static void _particleProblemSolveContacts(ParticleProblem *p, double dt, double
for (int i= 0; i < DIMENSION; ++i)
vloc[i] = p->velocity[c->o1 * DIMENSION + i] - p->segments[c->o0].v[i];
double pp = contactParticleBndSolve(c, vloc, vlocfree, &vn, dt);
//if (vn <= 0) pp = 0;
//else {
for (int i = 0; i < DIMENSION; ++i)
xn[i] = p->position[c->o1 * DIMENSION + i] + p->particles[c->o1].r * c->n[i] + vlocfree[i] * c->D / vn;
if (!segmentProjectionIsInside(&p->segments[c->o0], xn))
pp = 0;
//}
for (int i = 0; i < DIMENSION; ++i)
xn[i] = p->position[c->o1 * DIMENSION + i] + p->particles[c->o1].r * c->n[i] + vlocfree[i] * c->D / vn;
if (!segmentProjectionIsInside(&p->segments[c->o0], xn))
pp = 0;
double dp = pp - c->dv;
residual = fmax(fabs(dp), residual);
dp *= p->relax;
......@@ -522,13 +525,13 @@ double particleProblemMaxDt(const ParticleProblem *p, double alert) {
return alert / q;
}
void particleProblemIterate(ParticleProblem *p, double alert, double dt, double tol)
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 i = 0; i < DIMENSION; ++i)
p->velocity[j * DIMENSION + i] += p->externalForces[j * DIMENSION + i] * dt / p->particles[j].m;
_particleProblemGenContacts(p, alert);
_particleProblemSolveContacts(p, dt, tol);
_particleProblemSolveContacts(p, dt, tol, maxit);
for (size_t i = 0; i < vectorSize(p->position); ++i)
p->position[i] += p->velocity[i] * dt;
}
......@@ -668,6 +671,7 @@ int *particleProblemDiskTag(ParticleProblem *p) {return p->diskTag;}
int *particleProblemSegmentTag(ParticleProblem *p) {return p->segmentTag;};
#if DIMENSION == 3
int *particleProblemTriangleTag(ParticleProblem *p) {return p->triangleTag;};
double *particleProblemTriangle(ParticleProblem *p) {return (double*)&p->triangles[0];}
#endif
double *particleProblemVelocity(ParticleProblem *p) {return &p->velocity[0];}
double *particleProblemDisk(ParticleProblem *p) {return (double*)&p->disks[0];}
......
......@@ -8,7 +8,7 @@ ParticleProblem *particleProblemNew();
void particleProblemDelete(ParticleProblem *p);
void particleProblemLoad(ParticleProblem *p, const char *filename);
void particleProblemWrite(const ParticleProblem *p, const char *filename);
void particleProblemIterate(ParticleProblem *p, double alert, double dt, double tol);
void particleProblemIterate(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);
......
......@@ -80,7 +80,7 @@ struct ParticleProblem{};
void addBoundaryTriangle(const double x0[DIMENSION], const double x1[DIMENSION], const double x2[DIMENSION], int 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) { particleProblemIterate($self->_this, alert, dt, tol); }
void iterate(double alert, double dt, double tol = 0.005, int maxit = -1) { particleProblemIterate($self->_this, alert, dt, tol, maxit); }
double maxdt(double alert) { return particleProblemMaxDt($self->_this, alert); }
void load(const char *filename) { particleProblemLoad($self->_this, filename); }
......@@ -102,6 +102,12 @@ struct ParticleProblem{};
double *d = particleProblemSegment($self->_this);
return coordVector2PyArray(d, vectorSize(d) / (3 * DIMENSION), 3 * DIMENSION);
}
#if DIMENSION == 3
PyObject *triangles() {
double *d = particleProblemTriangle($self->_this);
return coordVector2PyArray(d, vectorSize(d) / (4 * DIMENSION), 4 * DIMENSION);
}
#endif
PyObject *position() {return coordVector2PyArray(particleProblemPosition($self->_this), particleProblemNParticle($self->_this), DIMENSION);}
PyObject *externalForces() {return coordVector2PyArray(particleProblemExternalForces($self->_this), particleProblemNParticle($self->_this), DIMENSION);}
void useJacobi(int jacobi, double relax = 0.5) {particleProblemUseJacobi($self->_this, jacobi, relax);}
......
......@@ -187,8 +187,8 @@ class ObjectCollection {
glClear(GL_COLOR_BUFFER_BIT|GL_DEPTH_BUFFER_BIT);
glColor4f(1, 0.5, 0, 1);
glPushMatrix();
double S = *std::max_element(radius.begin(), radius.end());
SolidSphere sphere(S, 20, 20);
double S = *std::max_element(radius.begin(), radius.end())/2;
SolidSphere sphere(S, 15, 15);
if (dim == 3) {
for(size_t i = 0; i < radius.size(); ++i) {
glPushMatrix();
......
from dgpy import *
import numpy as np
_U_ = 0
_V_ = 1
_W_ = 2
_UVW_=(0, 1, 2)
_UGT_ = (0, 4, 8)
_VGT_ = (1, 5, 9)
_WGT_ = (2, 6, 10)
_PGT_ = (3, 7, 11)
_P_ = 3
class brinkman(dgConservationLawPython) :
def _psiTerm0(self, cm, val, sol, solGrad, k, ca, vb):
k = k[:, np.newaxis]
ca = ca[:, np.newaxis]
invk = (1-ca)**2/(ca**2 *k)
cb = 1 - ca
invkb = (1-ca)/(ca * k)
dp = solGrad[:, 3 * _P_ : 3 * _P_ + 3]
du = solGrad[:, 3 * _U_ : 3 * _U_ + 3]
dv = solGrad[:, 3 * _V_ : 3 * _V_ + 3]
dw = solGrad[:, 3 * _W_ : 3 * _W_ + 3]
val[:, _UVW_] = - dp * ca - sol[:, _UVW_] * invk + (vb - self.stab * (1-ca) * dp) * invkb
val[:, _P_] = -du[:, 0] - dv[:, 1] - dw[:, 2]
def _psiTerm1(self, cm, val, sol, solGrad, mu, vb, ca):
mu = mu[:, np.newaxis]
ca = ca[:, np.newaxis]
u = sol[:, _U_]
v = sol[:, _V_]
w = sol[:, _W_]
du = solGrad[:, 3 * _U_ : 3 * _U_ + 3]
dv = solGrad[:, 3 * _V_ : 3 * _V_ + 3]
dw = solGrad[:, 3 * _W_ : 3 * _W_ + 3]
dp = solGrad[:, 3 * _P_ : 3 * _P_ + 3]
val[:, _UGT_] = -du * mu
val[:, _VGT_] = -dv * mu
val[:, _WGT_] = -dw * mu
s = self.stab * (1 - ca) + 1e-8
val[:, _PGT_] = - dp * s + vb
#if self.ns :
# val[:, _U_] += - u * u / ca
# val[:, _V_] += - u * v / ca
# val[:, _U_ + 3] += - v * u / ca
# val[:, _V_ + 3] += - v * v / ca
def _interfaceTerm(self, cm, val, sol, solGrad, n, mu, vb) :
ubn = vb[:, 0] * n[:, 0] + vb[:, 1] * n[:, 1] + vb[:, 2] * n[:, 2]
un = sol[:, _U_] * n[:, 0] + sol[:, _V_] * n[:, 1] + sol[:, W] * n[:, 2]
du = solGrad[:, 3 * _U_ : 3 * _U_ + 3]
dv = solGrad[:, 3 * _V_ : 3 * _V_ + 3]
dw = solGrad[:, 3 * _W_ : 3 * _W_ + 3]
dudn = du[:, 0] * n[:, 0] + du[:, 1] * n[:, 1] + du[:, 2] * n[:, 2]
dvdn = dv[:, 0] * n[:, 0] + dv[:, 1] * n[:, 1] + dv[:, 2] * n[:, 2]
dwdn = dw[:, 0] * n[:, 0] + dw[:, 1] * n[:, 1] + dw[:, 2] * n[:, 2]
val[:, _U_] = - mu * dudn
val[:, _V_] = - mu * dvdn
val[:, _W_] = - mu * dwdn
val[:, _P_] = 0 #-un + ubn * 0
def _velocity(self, cm, val, sol) :
val[:, 0] = sol[:, _U_]
val[:, 1] = sol[:, _V_]
val[:, 2] = sol[:, _W_]
def __init__(self, mu, k, porosity, particleMomentum, stab, ns) :
self.stab = stab
dgConservationLawPython.__init__(self, 4)
self.mu = mu
self.k = k
self.ns = ns
self.vb = particleMomentum
self.ca = porosity
self.t0 = functionNumpy(4, self._psiTerm0, [function.getSolution(), function.getSolutionGradient(), k, porosity, self.vb])
self.setVolumeTerm(0, self.t0)
self.t1 = functionNumpy(12, self._psiTerm1, [function.getSolution(), function.getSolutionGradient(), mu, self.vb, self.ca])
self.setVolumeTerm(1, self.t1)
self.b0 = functionNumpy(4, self._interfaceTerm, [function.getSolution(), function.getSolutionGradient(), function.getNormals(), mu, self.vb])
self.setInterfaceTerm(0, self.b0)
self.pressure = functionExtractCompNew(function.getSolution(), [_P_])
self.velocity = functionNumpy(3, self._velocity, [function.getSolution()])
self.setIsTemporal(_P_, False)
self.setIsLinear(not ns)
def setP1P1Mini(self, dof):
groups = dof.getGroups()
for i in range(groups.getNbElementGroups()) :
elements = groups.getElementGroup(i)
if elements.getFunctionSpace().type != MSH_TET_MINI :
Msg.Fatal("To use the p1 - p1 mini formulation, use p1 mini (i.e. order = -2) dgGroupCollection")
dof.removeDof(elements, _P_, 4)
#!/usr/bin/env python
import scontact3
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)
p = scontact3.ParticleProblem()
def genInitialPosition(p, N, rmin, rmax, L, H, rhop) :
n = 0
RM = rmax
R = rmin + (rmax - rmin ) * np.random.random([N])
for z in np.arange(H-RM, -H+RM, 2* -RM) :
for y in np.arange(L, -L, 2 * -RM) :
for x in np.arange(L, -L, 2 * -RM) :
R2 = x**2 + y**2
if R2 < (L - RM*1.5)**2:
p.addParticle((x, y, z), R[n], 4./3*R[n]**3 * np.pi * rhop);
n +=1
if n == N:
return
def loadMeshBoundary(p, filename, tags, shift = [0, 0, 0]):
import mesh
m = mesh.mesh(filename)
pnum = [m.getPhysicalNumber(2, i) for i in tags]
addv = set()
adde = set()
def vshift(v) :
return v[0] + shift[0], v[1] + shift[1], v[2] + shift[2]
for ent in m.entities :
if ent.dimension != 2 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(vshift(v), 0.000, tag)
for i in range(3) :
v0 = el.vertices[i]
v1 = el.vertices[(i + 1) % 3]
if set((v0[3], v1[3])) in adde :
continue
p.addBoundarySegment(vshift(v0), vshift(v1), tag)
p.addBoundaryTriangle(vshift(el.vertices[0]), vshift(el.vertices[1]), vshift(el.vertices[2]), tag)
r = 0.04
dr = 0.005
N = 8000
genInitialPosition(p, N = N, rmin = r -dr, rmax = r + dr , L = 1, H = 1.5, rhop=2)
p.useJacobi(False)
loadMeshBoundary(p, "mesh.msh", ["Box"])
g = 1
tEnd = 10
t = 0
iter = 0
tic = time.clock()
p.write("%s/part-%05i" % (outputdir, 0))
ioutput = 0
p.velocity()[:, :] = 0;
while t < tEnd :
p.externalForces()[:, :] = - p.velocity()
p.externalForces()[:, 2] = -g
p.externalForces()[:, :] *= p.particles()[:, [1]]
dt = p.maxdt(r +dr)/5
print(dt)
p.iterate(r + dr, dt, 1e-4)
t += dt
iter += 1
if iter %10 == 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, 0.5])
p.write("deposited")
class MshBaseElementType :
def __init__(self, name, tag, dimension) :
self.name = name
self.tag = tag
self.dimension = dimension
self.mshType = [None] * 11
self.mshTypeSerendip = [None] * 11
BaseType = [None]
Type = [None]
for t in [
('TYPE_PNT', 1, 0),
('TYPE_LIN', 2, 1),
('TYPE_TRI', 3, 2),
('TYPE_QUA', 4, 2),
('TYPE_TET', 5, 3),
('TYPE_PYR', 6, 3),
('TYPE_PRI', 7, 3),
('TYPE_HEX', 8, 3),
('TYPE_POLYG', 9, 2),
('TYPE_POLYH', 10, 3)
] :
BaseType.append(MshBaseElementType(*t))
globals()[t[0]] = BaseType[-1]
class MshElementType :
def __init__(self, name, tag, numVertices, baseType, order, serendip) :
self.name = name
self.tag = tag
self.numVertices = numVertices
self.baseType = baseType
self.order = order
self.serendip = serendip
for i, t in enumerate([
'MSH_LIN_2', 'MSH_TRI_3', 'MSH_QUA_4', 'MSH_TET_4', 'MSH_HEX_8',
'MSH_PRI_6', 'MSH_PYR_5', 'MSH_LIN_3', 'MSH_TRI_6', 'MSH_QUA_9',
'MSH_TET_10', 'MSH_HEX_27', 'MSH_PRI_18', 'MSH_PYR_14', 'MSH_PNT',
'MSH_QUA_8', 'MSH_HEX_20', 'MSH_PRI_15', 'MSH_PYR_13', 'MSH_TRI_9',
'MSH_TRI_10', 'MSH_TRI_12', 'MSH_TRI_15', 'MSH_TRI_15I', 'MSH_TRI_21',
'MSH_LIN_4', 'MSH_LIN_5', 'MSH_LIN_6', 'MSH_TET_20', 'MSH_TET_35',
'MSH_TET_56', 'MSH_TET_22', 'MSH_TET_28', 'MSH_POLYG_', 'MSH_POLYH_',
'MSH_QUA_16', 'MSH_QUA_25', 'MSH_QUA_36', 'MSH_QUA_12', 'MSH_QUA_16I',
'MSH_QUA_20', 'MSH_TRI_28', 'MSH_TRI_36', 'MSH_TRI_45', 'MSH_TRI_55',
'MSH_TRI_66', 'MSH_QUA_49', 'MSH_QUA_64', 'MSH_QUA_81', 'MSH_QUA_100',
'MSH_QUA_121', 'MSH_TRI_18', 'MSH_TRI_21I', 'MSH_TRI_24', 'MSH_TRI_27',
'MSH_TRI_30', 'MSH_QUA_24', 'MSH_QUA_28', 'MSH_QUA_32', 'MSH_QUA_36I',
'MSH_QUA_40', 'MSH_LIN_7', 'MSH_LIN_8', 'MSH_LIN_9', 'MSH_LIN_10',
'MSH_LIN_11', 'MSH_LIN_B', 'MSH_TRI_B', 'MSH_POLYG_B', 'MSH_LIN_C',
'MSH_TET_84', 'MSH_TET_120', 'MSH_TET_165', 'MSH_TET_220', 'MSH_TET_286',
None, None, None, 'MSH_TET_34', 'MSH_TET_40',
'MSH_TET_46', 'MSH_TET_52', 'MSH_TET_58', 'MSH_LIN_1', 'MSH_TRI_1',
'MSH_QUA_1', 'MSH_TET_1', 'MSH_HEX_1', 'MSH_PRI_1', 'MSH_PRI_40',
'MSH_PRI_75', 'MSH_HEX_64', 'MSH_HEX_125', 'MSH_HEX_216', 'MSH_HEX_343',
'MSH_HEX_512', 'MSH_HEX_729', 'MSH_HEX_1000', 'MSH_HEX_32', 'MSH_HEX_44',
'MSH_HEX_56', 'MSH_HEX_68', 'MSH_HEX_80', 'MSH_HEX_92', 'MSH_HEX_104',
'MSH_PRI_126', 'MSH_PRI_196', 'MSH_PRI_288', 'MSH_PRI_405', 'MSH_PRI_550',
'MSH_PRI_24', 'MSH_PRI_33', 'MSH_PRI_42', 'MSH_PRI_51', 'MSH_PRI_60',
'MSH_PRI_69', 'MSH_PRI_78', 'MSH_PYR_30', 'MSH_PYR_55', 'MSH_PYR_91',
'MSH_PYR_140', 'MSH_PYR_204', 'MSH_PYR_285', 'MSH_PYR_385', 'MSH_PYR_21',
'MSH_PYR_29', 'MSH_PYR_37', 'MSH_PYR_45', 'MSH_PYR_53', 'MSH_PYR_61',
'MSH_PYR_69', 'MSH_PYR_1', 'MSH_PNT_SUB', 'MSH_LIN_SUB', 'MSH_TRI_SUB',
'MSH_TET_SUB', 'MSH_TET_16'
]) :
if not t :
Type.append(None)
continue
w = t.split("_")
baseType = globals()["TYPE_"+w[1]]
nnodes2order = {
TYPE_LIN:(dict((o+1,o) for o in range(11)), {}),
TYPE_TRI:(dict(((o+1) * (o+2) / 2,o) for o in range(11)), dict((3 * o,o) for o in range(11))),
TYPE_QUA:(dict(((o+1) * (o+1),o) for o in range(11)), dict((4 * o,o) for o in range(11))),
TYPE_TET:(dict(((o+1) * (o+2) * (o+3) / 6,o) for o in range(11)), dict((4 + 6 * (o - 1), o) for o in range(11))),
TYPE_HEX:(dict(((o+1) * (o+1) * (o+1),o) for o in range(11)), dict((8 + 12 * (o - 1), o) for o in range(11))),
TYPE_PRI:(dict(((o+1) * (o+2) * (o+1)/2,o) for o in range(11)), dict((6 + 9 * (o - 1), o) for o in range(11))),
TYPE_PYR:(dict(((o+1) * (o+2) * (2*o +3)/6, o) for o in range(11)), dict((5 + 8 * (o - 1),o) for o in range(11)))
}
try :
numnodes = int(w[2].strip("I"))
except:
numnodes = 1 if baseType == TYPE_PNT else -1
try :
serendip = (w[2][-1] == 'I')
except :
serendip = False
if numnodes == -1 :
order = -1
elif baseType == TYPE_PNT :
order = 0
serendip = False
else:
if (not serendip) and numnodes in nnodes2order[baseType][0] :
order = nnodes2order[baseType][0][numnodes]
else :
order = nnodes2order[baseType][1][numnodes]
serendip = True
Type.append(MshElementType(t, i + 1, numnodes, baseType, order, serendip))
if serendip :
baseType.mshTypeSerendip[order] = Type[-1]
else :
baseType.mshType[order] = Type[-1]
globals()[t] = Type[-1]
L = 1;
l = 0.1;
H = 1.5;
Lc = 0.3;
h = 0.5;
lc = 0.1;
Point(1) = {-L, 0, H, Lc};
Point(2) = {0, -L, H, Lc};
Point(3) = {L, 0, H, Lc};
Point(4) = {0, L, H, Lc};
Point(5) = {0, 0, H, Lc};
Circle(1) = {1, 5, 2};
Circle(2) = {2, 5, 3};
Circle(3) = {3, 5, 4};
Circle(4) = {4, 5, 1};
Line Loop(1) = {1, 2, 3, 4};
Plane Surface(1) = {1};
Point(6) = {-l, 0, h, lc};
Point(7) = {0, -l, h, lc};
Point(8) = {l, 0, h, lc};
Point(9) = {0, l, h, lc};
Point(10) = {0, 0, h, lc};
Circle(6) = {6, 10, 7};
Circle(7) = {7, 10, 8};
Circle(8) = {8, 10, 9};
Circle(9) = {9, 10, 6};
Line Loop(2) = {6, 7, 8, 9};
Plane Surface(2) = {2};
tmp[] = Extrude {0, 0, -2*H}{
Surface{1};
};
tmp2[] = Extrude {0, 0, -2*h}{
Surface{2};
};
Delete {Volume {tmp[1],tmp2[1]};}
Surface Loop(10) = {1, tmp[0], tmp[2], tmp[3], tmp[4], tmp[5]};
Surface Loop(11) = {2, tmp2[0], tmp2[2], tmp2[3], tmp2[4], tmp2[5]};
Volume(10) = {10, 11};
Physical Volume("Domain") = {10};
Physical Surface("Box") = {1, tmp[0], tmp[2], tmp[3], tmp[4], tmp[5]};
Physical Surface("Pile") = {2, tmp2[0], tmp2[2], tmp2[3], tmp2[4], tmp2[5]};
Physical Point("PtFix") = {1};
import gmshType
import time
class mesh() :
class element():
def __init__(self, etype, vertices, partition, tag):
self.etype = etype
self.vertices = vertices
self.partition = partition
self.tag = tag
def genFaces(self) :
v = self.vertices
order = self.etype.order
if self.etype.baseType == gmshType.TYPE_TRI :
return [mesh.element(gmshType.TYPE_LIN.mshType[order], [v[i], v[(i + 1) % 3]] + v[3 + i * (order-1):3 + (i+1) * (order-1)], self.partition, -1) for i in range(3)]
elif self.etype.baseType == gmshType.TYPE_QUA :
return [mesh.element(gmshType.TYPE_LIN.mshType[order], [v[i], v[(i + 1) % 4]] + v[4 + i * (order-1):4 + (i+1) * (order-1)], self.partition, -1) for i in range(4)]
elif self.etype.baseType == gmshType.TYPE_LIN :
return [mesh.element(gmshType.MSH_PNT, [v[i]], self.partition, -1) for i in range(2)]
class entity():
def __init__(self, tag, dim, physicals) :
self.tag = tag
self.dimension = dim
self.physicals = physicals
self.elements = []
def __init__(self, filename = None) :
self.vertices = []
self.physicals = [{}, {}, {}, {}]
self.entities = []
self.maxeid = 0
self.useFormat3 = False
if filename :
self.read(filename)
def write(self, filename, format3=None) :
output = open(filename, "w")
if format3 == None :
format3 = self.useFormat3
if format3 :
output.write("$MeshFormat\n3 0 8\n$EndMeshFormat\n")
output.write("$PhysicalNames\n%i\n" % (len(self.physicals[0]) + len(self.physicals[1]) + len(self.physicals[2]) + len(self.physicals[3])))
for dim, plist in enumerate(self.physicals) :
for name, tag in plist.items() :
output.write("%i %i \"%s\"\n" % (dim, tag, name))
output.write("$EndPhysicalNames\n")
output.write("$Entities\n%i\n" % (len(self.entities)))
for e in self.entities :
output.write("%i %i %i %s\n" % (e.tag, e.dimension, len(e.physicals), " ".join(repr(ip) for ip in e.physicals)))
output.write("$EndEntities\n")
output.write("$Nodes\n%d\n" % len(self.vertices))
for (x, y, z, j) in self.vertices :
output.write("%i %.16g %.16g %.16g 0\n" % (j, x, y, z))
output.write("$EndNodes\n")
output.write("$Elements\n%d\n" % sum ([len(e.elements) for e in self.entities]))
for entity in self.entities :
for e in entity.elements :
ntag = len(e.vertices) + ((1 + len(e.partition)) if e.partition else 0)
spart = " ".join(repr(i) for i in [len(e.partition)] + e.partition) if e.partition else ""
output.write("%i %i %i %i %s %s\n" % (e.tag, e.etype.tag, entity.tag, ntag, " ".join(repr(v[3]) for v in e.vertices), spart))
output.write("$EndElements\n")
else :
output.write("$MeshFormat\n2.2 0 8\n$EndMeshFormat\n")
output.write("$PhysicalNames\n%i\n" % (len(self.physicals[0]) + len(self.physicals[1]) + len(self.physicals[2]) + len(self.physicals[3])))
for dim, plist in enumerate(self.physicals) :
for name, tag in plist.items() :
output.write("%i %i \"%s\"\n" % (dim, tag, name))
output.write("$EndPhysicalNames\n")
output.write("$Nodes\n%d\n" % len(self.vertices))
for (x, y, z, j) in self.vertices :
output.write("%i %.16g %.16g %.16g\n" % (j, x, y, z))
output.write("$EndNodes\n")
output.write("$Elements\n%d\n" % sum ([len(e.elements) for e in self.entities]))
for entity in self.entities :
for e in entity.elements :
ntag = 2 + ((1 + len(e.partition)) if e.partition else 0)
spart = " ".join(repr(i) for i in [len(e.partition)] + e.partition) if e.partition else ""
ptag = entity.physicals[0] if entity.physicals else 0
output.write("%i %i %i %i %i %s %s\n" % (e.tag, e.etype.tag, ntag, ptag, entity.tag, spart, " ".join(repr(v[3]) for v in e.vertices)))
output.write("$EndElements\n")
def read(self, filename):
fin = open(filename, "r");
l = fin.readline()
vmap = {}
entitymap = {}
while l != "" :
w = l.split()
if w[0] == "$MeshFormat":
l = fin.readline().split()
if float(l[0]) == 3.:
self.useFormat3 = True
elif int(float(l[0])) == 2 :
self.useFormat3 = False
else :
print("error : cannot read mesh format " + l[0])
l = fin.readline()
elif w[0] == "$PhysicalNames" :
n = int(fin.readline())
for i in range(n) :
dim, tag, name = fin.readline().split()
self.physicals[int(dim)][name[1:-1]] = int(tag)
fin.readline()
elif w[0] == "$Entities" and self.useFormat3:
n = int(fin.readline())
for i in range(n) :
l = fin.readline().split()
j, dim, nphys = int(l[0]), int(l[1]), int(l[2])
self.entities.append(mesh.entity(j, dim, [int(ip) for ip in l[3:3+nphys]]))
entitymap[(dim, j)] = self.entities[-1]
fin.readline()
elif w[0] == "$Nodes" :
n = int(fin.readline())
for i in range(n) :
if self.useFormat3 :
(j, x, y, z, t) = fin.readline().split()
else :
(j, x, y, z) = fin.readline().split()
self.vertices.append([float(x), float(y), float(z), int(j)])
vmap[int(j)] = self.vertices[-1]
elif w[0] == "$Elements" :
n = int(fin.readline())
for i in range(n) :
l = fin.readline().split()
if self.useFormat3 :
j, t, e, nf = int(l[0]), int(l[1]), int(l[2]), int(l[3])
nv = gmshType.Type[t].numVertices
vertices = [vmap[int(i)] for i in l[4:4+nv]]
partition = [int(i) for i in l[5 + nv : 5 + nv + int(l[4 + nv])]] if nf > nv else []
else :
j, t, nf, p, e = int(l[0]), int(l[1]), int(l[2]), int(l[3]), int(l[4])
vertices = [vmap[int(i)] for i in l[3 + nf:]]
partition = [int(i) for i in l[6 : 6 + int(l[5])]] if nf > 2 else []
edim = gmshType.Type[t].baseType.dimension
entity = entitymap.get((edim, e), None)
if not self.useFormat3 and not entity:
entity = mesh.entity(e, edim, [p])
self.entities.append(entity)
entitymap[(edim, e)] = entity
entity.elements.append(mesh.element(gmshType.Type[t], vertices, partition, j))
self.maxeid = max(self.maxeid, j)
l = fin.readline()
def getPhysicalNumber(self, dim, name) :
t = self.physicals[dim].get(name, None)
if not t :
self.physicals[dim][name] = max(self.physicals[dim].values()) + 1 if self.physicals[dim] else 1
t = self.physicals[dim][name]
return t