How can I create a triangle mesh file for ESyS-Particle simulations?

ESyS-Particle provides a facility to define walls or obstacles of non-planar geometry via triangular meshes. In this way complex geometries such as cylindrical containers for particles can be simulated. Triangular meshes are imported into ESyS-Particle simulations via LsmMpi.readMesh(..) and need to be provided as an ASCII text file employing a bespoke mesh file format. The ESyS-Particle Tutorial (https://launchpad.net/esys-particle/2.3/2.3.1/+download/ESyS-Particle_Tutorial.pdf) contains basic instructions on how to construct a suitable mesh file and import it into a simulation.

Whilst it is relatively easy to define simple meshes by hand, more complex meshes are best constructed using CAD software. A variety of both proprietary and open source CAD packages may be used. ESyS-Particle developers use and recommend OpenSCAD (https://www.openscad.org/). Regardless of the CAD package used, the mesh geometry should be exported in Stereolithograph (STL) format. Thereafter, use the following script (TriMeshTools.py) to convert the STL mesh into ESyS-Particle's mesh file format.

You will need to install the Python bindings to the VTK library before using this script. To do so on Ubuntu-20.04 or later, execute the following command from the terminal:

sudo apt install python3-vtk

#BEGIN: TriMeshTools.py

from gengeo import *
from math import *
from random import randint
from vtk import *

def readTriMeshFile(inputfile="Sample_0.msh"):
   infile = open(inputfile,"r")
   lines = infile.readlines()
   infile.close()

   npts = int(lines[1].split(' ')[1])
   pts = []
   Xmin = 1.e10
   Xmax = -1.e10
   Ymin = 1.e10
   Ymax = -1.e10
   Zmin = 1.e10
   Zmax = -1.e10
   for line in lines[2:2+npts]:
      data = line.split(' ')
      X = float(data[3])
      Y = float(data[4])
      Z = float(data[5])
      if (X<Xmin): Xmin = X
      if (X>Xmax): Xmax = X
      if (Y<Ymin): Ymin = Y
      if (Y>Ymax): Ymax = Y
      if (Z<Zmin): Zmin = Z
      if (Z>Zmax): Zmax = Z
      pts.append(Vector3(X,Y,Z))

   ntris = int(lines[3+npts].split(' ')[1])
   tris = []
   for line in lines[4+npts:]:
      data = line.split(' ')
      p1 = int(data[2])
      p2 = int(data[3])
      p3 = int(data[4])
      tris.append([p1,p2,p3,2])

   bbox = [Vector3(Xmin,Ymin,Zmin),Vector3(Xmax,Ymax,Zmax)]

   return pts,tris,bbox

def writeTriMeshFile (pts,tris,outputfile="Sample_0.msh"):
   outstr = "Triangle\n3D-Nodes %d\n" % (len(pts))
   pid = 0
   for pt in pts:
      X = pt.X()
      Y = pt.Y()
      Z = pt.Z()
      outstr += "%d %d 0 %.10e %.10e %.10e\n" % (pid,pid,X,Y,Z)
      pid += 1

   outstr += "\nTri3 %d\n" % (len(tris))
   tid = 0
   for tri in tris:
      outstr += "%d 0 %d %d %d\n" % (tid,tri[0],tri[1],tri[2])

   outfile = open(outputfile,"w")
   outfile.write(outstr)
   outfile.close()
   return

def translate (pts,tris,bbox,displacement):
   dpts = []
   for pt in pts:
      dpts.append(pt+displacement)

   dbbox = [bbox[0]+displacement,bbox[1]+displacement]

   return dpts,tris,dbbox

def scale (pts,tris,bbox,scaleFactor=0.5,newCentrePt=None):
   spts = []
   minPt = bbox[0]
   maxPt = bbox[1]
   centrePt = (maxPt-minPt)*0.5 + minPt
   dpts,dtris,dbbox = translate(pts,tris,bbox,centrePt*-1)
   spts = []
   for pt in dpts:
      spts.append(pt*scaleFactor)
   sbbox = [dbbox[0]*scaleFactor,dbbox[1]*scaleFactor]
   if (newCentrePt == None):
      opts,otris,obbox = translate(spts,tris,sbbox,centrePt)
   else:
      opts,otris,obbox = translate(spts,tris,sbbox,newCentrePt)

   return opts,otris,obbox

def scale_and_translate (pts,tris,bbox,maxLength=10.0,centrePt=Vector3(0,0,0)):
   dbox = bbox[1] - bbox[0]
   Lx = dbox.X()
   Ly = dbox.Y()
   Lz = dbox.Z()
   currMaxLength = max([Lx,Ly,Lz])
   scaleFactor = maxLength/currMaxLength

   newpts,newtris,newbbox = scale(pts,tris,bbox,scaleFactor=scaleFactor,newCentrePt=centrePt)
   return newpts,newtris,newbbox

def normalVector(p1,p2,p3):
   r1 = [p2[0]-p1[0],p2[1]-p1[1],p2[2]-p1[2]]
   r2 = [p3[0]-p1[0],p3[1]-p1[1],p3[2]-p1[2]]

   cp = [r2[1]*r1[2] - r2[2]*r1[1], r2[2]*r1[0] - r2[0]*r1[2], r2[0]*r1[1] - r2[1]*r1[0]]

   norm = sqrt(cp[0]*cp[0]+cp[1]*cp[1]+cp[2]*cp[2])

   if (norm > 0.0):
      nV = [cp[0]/norm,cp[1]/norm,cp[2]/norm]
   else:
      nV = [0,0,0]
      print ("yikes!")

   return nV

def reorderNodes(pts, tris, inPoint = [0.,0.,0.3]):
   for tid in range(len(tris)):
      id1 = tris[tid][0]
      id2 = tris[tid][1]
      id3 = tris[tid][2]
      p1 = [pts[id1].X(), pts[id1].Y(), pts[id1].Z()]
      p2 = [pts[id2].X(), pts[id2].Y(), pts[id2].Z()]
      p3 = [pts[id3].X(), pts[id3].Y(), pts[id3].Z()]

      nV = normalVector(p1,p2,p3)

      dp = [inPoint[0]-p1[0],inPoint[1]-p1[1],inPoint[2]-p1[2]]

      dp_nV = dp[0]*nV[0] + dp[1]*nV[1] + dp[2]*nV[2]

      if (dp_nV < 0):
         tris[tid][0] = id3
         tris[tid][2] = id1

   return

def toTriPatchSet(pts,tris,bbox):
      tripatchset = TriPatchSet()
      for tri in tris:
         p1 = pts[tri[0]]
         p2 = pts[tri[1]]
         p3 = pts[tri[2]]
         tag = tri[3]
         tripatchset.addTriangle(p1,p2,p3,tag)
      return tripatchset,bbox

def readstl(inputfile="mill_liner.stl"):
   reader = vtkSTLReader()
   reader.SetFileName (inputfile)
   reader.Update()

   thePoly = vtkPolyData()
   thePoly.ShallowCopy(reader.GetOutput())

   pts = []
   for i in range (thePoly.GetNumberOfPoints()):
      p = thePoly.GetPoint(i)
      pts.append(Vector3(p[0],p[1],p[2]))

   tris = []
   for i in range (thePoly.GetNumberOfPolys()):
      t = thePoly.GetCell(i)
      id0 = t.GetPointIds().GetId(0)
      id1 = t.GetPointIds().GetId(1)
      id2 = t.GetPointIds().GetId(2)
      tris.append([id0,id1,id2,2])

   bounds = thePoly.GetBounds()
   minPt = Vector3(bounds[0],bounds[2],bounds[4])
   maxPt = Vector3(bounds[1],bounds[3],bounds[5])
   bbox = [minPt,maxPt]

   return pts,tris,bbox

if __name__ == "__main__":
   pts,tris,bbox = readstl("test.stl")
   pts,tris,bbox = scale_and_translate(pts,tris,bbox,maxLength=18.0,centrePt=Vector3(0,0,0))
   writeTriMeshFile(pts,tris,"test.msh")

#END: TriMeshTools.py

Older Notes (now Deprecated):
=========================

To construct suitable mesh walls for ESyS-Particle simulations, the developers recommend the third party, Open Source mesh generation software called gmsh:
http://gmsh.info/

To learn how to use gmsh to construct a simple cylindrical mesh, the following step-by-step guide may be helpful (albeit illustrating usage of an older version of gmsh):
http://www.esys-particle.org/files/Gmsh_tutorial_for_ESyS-Particle.pdf
There are also other gmsh tutorials available online.

For gmsh version 3.0 and later, save your triangle mesh in gmsh (.msh) format then use the following Python script to convert the file into the ESyS-Particle trimesh file format:
http://www.esys-particle.org/wiki/pmwiki.php/Main/Gmsh2TriMesh

[For older versions of gmsh, the following script will convert from 'Medit Mesh' format to ESyS-Particle trimesh format. This method is now deprecated]
http://www.esys-particle.org/wiki/pmwiki.php/Main/Medit2TriMesh

The TriMesh file created (cylinder.lsm) can then be imported into ESyS-Particle simulations using the LsmMpi.readMesh(..) subroutine. For more instructions on importing TriMeshes, please refer to the ESyS-Particle Tutorial:
https://launchpad.net/esys-particle/2.3/2.3.1/+download/ESyS-Particle_Tutorial.pdf

This presentation also explains how to import and use triangle meshes in ESyS-Particle:
http://www.esys-particle.org/files/TriangleMeshesOverview.pdf

N.B. The conversion scripts (Gmsh2TriMesh.py or medit2TriMesh.py) may require alteration for converting meshes with geometries differing from that in the step-by-step guide. Probably the most important change will be to specify the coordinates of a point lying "inside" the mesh (called inPoint in the scripts). Since ESyS-Particle TriMeshes are one-sided, this point needs to be specified on the same side of the mesh as the intended location of particles during simulations.