Commit 4764feb4 authored by Jonathan Lambrechts's avatar Jonathan Lambrechts
Browse files

use gmsh api for load_msh

parent 17001cc8
Pipeline #8681 failed with stages
in 1 minute and 18 seconds
......@@ -517,4 +517,9 @@ void mesh_set_elements(Mesh *m, int n_nodes, double *x, int n_elements, int *ele
for (int i = 0; i < n_physicals; ++i) {
m->boundary_names[i] = strdup(physicals[i]);
}
m->periodic_mapping = malloc(sizeof(int)*m->n_nodes);
m->n_periodic = 0;
for(int i = 0; i < m->n_nodes; ++i){
m->periodic_mapping[i] = i;
}
}
......@@ -33,6 +33,7 @@ import numpy as np
import signal
import os
import sys
import gmsh
from . import VTK
try :
from .petsclsys import LinearSystem
......@@ -43,6 +44,8 @@ except :
dir_path = os.path.dirname(os.path.realpath(__file__))
lib2 = np.ctypeslib.load_library("libmbfluid2",dir_path)
lib3 = np.ctypeslib.load_library("libmbfluid3",dir_path)
gmsh.initialize()
gmsh.option.setNumber("General.Terminal",1)
signal.signal(signal.SIGINT, signal.SIG_DFL)
......@@ -127,7 +130,45 @@ class FluidProblem :
Keyword argument:
mesh_file_name -- Name of the mesh.msh file containing information about the domain
"""
self._lib.fluid_problem_load_msh(self._fp, mesh_file_name.encode())
dim = self.dimension()
etype = 2 if dim == 2 else 4
betype = 1 if dim == 2 else 2
gmsh.open(mesh_file_name)
gmsh.model.mesh.renumberNodes()
_, x, _ = gmsh.model.mesh.getNodes()
x = x.reshape([-1,3])
_, el = gmsh.model.mesh.getElementsByType(etype)
el = (el-1).reshape([-1,dim+1])
bnd = []
btag = []
bname = []
for ip,(pdim, ptag) in enumerate(gmsh.model.getPhysicalGroups(dim-1)) :
bname.append(gmsh.model.getPhysicalName(pdim,ptag))
for etag in gmsh.model.getEntitiesForPhysicalGroup(pdim, ptag):
for eltype, _, elnodes in zip(*gmsh.model.mesh.getElements(pdim,etag)):
if eltype != betype:
continue
bnd.append(elnodes-1)
btag.append(np.full((elnodes.shape[0]//2), ip))
bnd = np.hstack(bnd).reshape([-1,dim])
btag = np.hstack(btag)
cbname = (c_char_p*len(bname))(*(i.encode() for i in bname))
#remove nodes and boundaries not connected to elements
keepnodes = np.full(x.shape[0], False, np.bool)
keepnodes[el] = True
newidx = np.cumsum(keepnodes)-1
el = newidx[el]
x = x[keepnodes,:]
keepbnd = np.full(bnd.shape[0], True, np.bool)
for i in range(dim) :
keepbnd = np.logical_and(keepbnd, keepnodes[bnd[:,i]])
bnd = newidx[bnd][keepbnd,:]
btag = btag[keepbnd]
self._lib.fluid_problem_set_elements(self._fp,
c_int(x.shape[0]),_np2c(x,np.float64),
c_int(el.shape[0]),_np2c(el,np.int32),
c_int(bnd.shape[0]),_np2c(bnd,np.int32),
_np2c(btag,np.int32),c_int(len(cbname)),cbname)
self.sys = None
def set_wall_boundary(self, tag, pressure=None, velocity=None) :
......
......@@ -116,7 +116,7 @@ fluid.set_wall_boundary("Top", pressure=0)
fluid.set_wall_boundary("Bottom")
fluid.set_wall_boundary("Lateral")
# Set location of the particles in the mesh and compute the porosity in each computation cell
fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity(),p.contact_forces(),init=True)
fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity(),p.contact_forces())
t = 0
fluid.solution()[:,2] = (fluid.coordinates()[:,1]-0.6)*rho*g[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