Commit 40c184d9 authored by Michel Henry's avatar Michel Henry
Browse files

Periodique 3D

parent a4e056b7
......@@ -308,6 +308,7 @@ static void _gmsh_mesh_read_4(GmshMesh *mesh, FileReader *r) {
int n_nodes_by_element;
if (element_type == 1) n_nodes_by_element = 2;
else if (element_type == 2) n_nodes_by_element = 3;
else if (element_type == 4) n_nodes_by_element = 4;
else if (element_type == 15) n_nodes_by_element = 1;
else file_reader_error_at(r,"unknown element type '%i'",element_type);
size_t n_element_block = file_reader_get_unsigned_int(r);
......@@ -318,7 +319,7 @@ static void _gmsh_mesh_read_4(GmshMesh *mesh, FileReader *r) {
}
}
if (entity == NULL) {
file_reader_error_at(r,"unkonwn entity");
file_reader_error_at(r,"unknown entity");
}
if (entity->n_nodes_by_element != 0 && entity->n_nodes_by_element != n_nodes_by_element) {
file_reader_error_at(r,"mixed meshes not suppoted");
......@@ -356,9 +357,9 @@ static void _gmsh_mesh_read_4(GmshMesh *mesh, FileReader *r) {
file_reader_next_line(r);
p->n_affine = file_reader_get_int(r);
p->affine = malloc(sizeof(double)*p->n_affine);
file_reader_next_line(r);
for (int j = 0; j < p->n_affine; ++j)
p->affine[j] = file_reader_get_double(r);
file_reader_next_line(r);
p->n_nodes = file_reader_get_int(r);
file_reader_next_line(r);
p->nodes = malloc(sizeof(size_t)*2*p->n_nodes);
......@@ -786,6 +787,7 @@ size_t gmsh_mesh_n_elements(GmshMesh* mesh, const int dim) {
for(int ie = 0; ie < mesh->n_entities[dim]; ++ie){
const Entity *e = mesh->entities[dim][ie];
if (e->n_nodes_by_element != dim +1) {
printf("Unvalid number by element\n");
exit(1);
}
n_elements += e->n_elements;
......
......@@ -442,7 +442,7 @@ class ParticleProblem :
def add_boundary_periodic_triangle(self, x0, x1, x2, etag):
if self._dim != 3 :
raise NameError("Triangle boundaries only available in 3D")
selb._lib.particleProblemAddBoundaryPeriodicTriangle.restype = c_size_t
self._lib.particleProblemAddBoundaryPeriodicTriangle.restype = c_size_t
return self._lib.particleProblemAddBoundaryPeriodicTriangle(self._p, self._coord_type(*x0), self._coord_type(*x1), self._coord_type(*x2), c_int(etag))
......@@ -615,7 +615,7 @@ class ParticleProblem :
if num == number :
return name
return ""
def get_entities(self,_mesh, dim=2):
def get_entities(self,_mesh, dim):
entities = []
_nnodes = self.n_nodes(_mesh, dim)
_nodes = self.nodes(_mesh, dim, _nnodes)
......
......@@ -119,6 +119,6 @@ while t < tEnd :
if ii %outf == 0 :
ioutput = int(ii/outf) + 1
p.write_vtk(outputdir, ioutput, t)
#fluid.export_vtk(outputdir, t, ioutput)
fluid.export_vtk(outputdir, t, ioutput)
ii += 1
print("%i : %.2g/%.2g (cpu %.6g)" % (ii, t, tEnd, time.time() - tic))
......@@ -15,4 +15,4 @@ Physical Surface("Z") = {5,6};
Physical Point("PtFix") = {5};
Characteristic Length{ PointsOf{ Volume{:}; } } = lc;
Mesh.MshFileVersion = 2;
//Mesh.MshFileVersion = 2;
......@@ -66,14 +66,16 @@ if not os.path.isdir(outputdir) :
os.makedirs(outputdir)
# Physical parameters
g = np.array([0,-9.81,0]) # gravity
# g = np.array([0,-9.81,0]) # gravity
g = np.array([0,0,0])
rho = 1000 # fluid density
rhop = 2640 # particles density
nu = 1e-6 # kinematic viscosity
L = 1
# Numerical parameters
outf = 50 # number of iterations between output files
dt = 1e-4 # time step
outf = 5 # number of iterations between output files
dt = 0.1 # time step
tEnd = 100 # final time
#
......@@ -81,10 +83,15 @@ tEnd = 100 # final time
#
# Particles structure creation and setting particles properties
r = 1e-2
p = scontact.ParticleProblem(3)
p.load_msh_boundaries("mesh.msh", ["Bottom", "Top", "X", "Z"])
genInitialPosition(p,outputdir, rhop)
print(p.n_particles())
p.add_particle((r,0.5,0.5), r, 4/3*r**3*np.pi*rhop)
p.add_particle((0.5,0.5,0.5), r, 4/3*r**3*np.pi*rhop)
state = p.state()
state.v[1, 0] += 0.01
p.set_state(state)
# genInitialPosition(p,outputdir, rhop)
p.write_vtk(outputdir,0,0)
# Initial time and iteration
......@@ -113,12 +120,19 @@ G[:,1] = p.mass()[:,0]*g[1]
# COMPUTATION LOOP
#
while t < tEnd :
time_integration.iterate(fluid, p, dt,min_nsub=20, contact_tol = 1e-8,external_particles_forces=G)
time_integration.iterate(None, p, dt,min_nsub=20, contact_tol = 1e-8,external_particles_forces=G)
state = p.state()
state.x[:,0] = np.fmod(state.x[:,0],L)
ind = np.where(state.x[:,0] <= 0)
state.x[ind,0] += L
p.set_state(state)
t += dt
# Output files writing
if ii %outf == 0 :
ioutput = int(ii/outf) + 1
p.write_vtk(outputdir, ioutput, t)
#fluid.export_vtk(outputdir, t, ioutput)
fluid.export_vtk(outputdir, t, ioutput)
ii += 1
print("%i : %.2g/%.2g (cpu %.6g)" % (ii, t, tEnd, time.time() - tic))
SetFactory("OpenCASCADE");
L = 0.014;
H = 0.036;
P = 0.01675;
L = 1;
H = 1;
P = 1;
y = 0;
lc = 1e-2;//2e-3;
lc= 1e-1;
SetFactory("OpenCASCADE");
Box(1) = {-L,-H,-P,2*L,2*H,2*P};
Periodic Surface {1} = {2};
Volume("Domain") = {1};
Box(1) = {0,0,0,L,H,P};
Periodic Surface {1} = {2} Translate {-L,0,0};
Physical Volume("Domain") = {1};
Physical Surface("Top") = {4};
Physical Surface("Bottom") = {3};
Physical Surface("X") = {1,2};
Physical Surface("Z") = {5,6};
Physical Point("PtFix") = {5};
Characteristic Length{ PointsOf{ Volume{:}; } } = lc;
Mesh.MshFileVersion = 2;
//Mesh.MshFileVersion = 2;
......@@ -40,7 +40,7 @@ outputdir = "outputDepotNP"
if not os.path.isdir(outputdir) :
os.makedirs(outputdir)
subprocess.call(["gmsh", "-2", "mesh.geo","-clscale","1"])
subprocess.call(["gmsh", "-2", "meshDepot.geo","-clscale","1"])
t = 0
ii = 0
......@@ -66,10 +66,10 @@ alpha = 1e-4 # stabilization coefficient
epsilon = alpha*lcmin**2 /nu # stabilization parametre
print('epsilon',epsilon)
shutil.copy("mesh.msh", outputdir +"/mesh.msh")
shutil.copy("meshDepot.msh", outputdir +"/mesh.msh")
p = scontact.ParticleProblem(2)
p.load_msh_boundaries("mesh.msh", ["Bottom", "Top", "Right", "Left"])
p.load_msh_boundaries("meshDepot.msh", ["Bottom", "Top", "Right", "Left"])
for x in np.arange(2*r,L-2*r,2.1*r):
for y in np.arange(H-0.2,H-r,2.1*r):
R = r*(1-np.random.random()/4)
......@@ -83,7 +83,7 @@ outf = 5 # number of iterations between ou
#Object fluid creation + Boundary condition of the fluid (field 0 is horizontal velocity; field 1 is vertical velocity; field 2 is pressure)
fluid = mbfluid.FluidProblem(2,g,nu*rho,rho,petsc_solver_type="-pc_type lu")
fluid.load_msh("mesh.msh")
fluid.load_msh("meshDepot.msh")
#fluid.set_open_boundary("Right",pressure = 0)
#fluid.set_open_boundary("Left", pressure = 0)
fluid.set_wall_boundary("Bottom",velocity=[0,0])
......
L = 0.4;
H = 0.6;
L = 1;
H = 1;
y = 0;
lc = 0.01;
lc = 0.1;
Point(1) = {0,0,0,lc};
Point(2) = {L,0,0,lc};
......
L = 0.4;
H = 0.6;
y = 0;
lc = 0.01;
Point(1) = {0,0,0,lc};
Point(2) = {L,0,0,lc};
Point(3) = {L,H,0,lc};
Point(4) = {0,H,0,lc};
Line(1) = {1,2};
Line(2) = {2,3};
Line(3) = {4,3};
Line(4) = {1,4};
Line Loop(1) = {1,2,-3,-4};
Plane Surface(1) = {1};
Periodic Curve {2} = {4};
Physical Line("Left") = {4};
Physical Line("Right") = {2};
Physical Line("Bottom") = {1};
Physical Line("Top") = {3};
Physical Surface("Domain") = {1};
Physical Point("PtFix") = {1};
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