#############################################################
##							   ##
## Copyright (c) 2003-2009 by The University of Queensland ##
## Earth Systems Science Computational Centre (ESSCC)	   ##
## http://www.uq.edu.au/esscc				   ##
##							   ##
## Primary Business: Brisbane, Queensland, Australia	   ##
## Licensed under the Open Software License version 3.0	   ##
## http://www.opensource.org/licenses/osl-3.0.php	   ##
##							   ##
#############################################################

# medit2TriMesh.py: converts Medit meshes to ESyS-Particle TriMesh format
#
#	Author: D. Weatherley
#	Organisation: Earth Systems Science Computational Centre
#		      The University of Queensland, Brisbane, Australia.
#	Date: 18 September 2009


#import required modules:
from math import sqrt
from sys import argv

def readMeditMesh (meshfile="corner.mesh"):
   infile = open(meshfile,"r")
   lines = infile.readlines()
   infile.close()

   nodelist = []
   numNodes = int(lines[4].split()[0])
   for i in range(numNodes):
      nodelist.append([0,0.0,0.0,0.0])

   id = 0
   for line in lines[5:numNodes+5]:
      data = line.split()
      X = float(data[0])
      Y = float(data[1])
      Z = float(data[2])
      nodelist[id][0] = id
      nodelist[id][1] = X
      nodelist[id][2] = Y
      nodelist[id][3] = Z
      id += 1

   numTriangles = int(lines[numNodes+6].split()[0])
   triList = []
   for i in range(numTriangles):
      triList.append([0,0,0,0])

   tid = 0
   for line in lines[numNodes+7:numNodes+7+numTriangles]:
      data = line.split()
      id1 = int(data[0])-1
      id2 = int(data[1])-1
      id3 = int(data[2])-1
      triList[tid][0] = tid
      triList[tid][1] = id1
      triList[tid][2] = id2
      triList[tid][3] = id3
      tid += 1

   return (nodelist, triList)

def writeTriMesh (nodeList, triList, outfilename="trimesh.lsm"):
   outstring = "Triangle\n3D-Nodes %d\n" % (len(nodeList))
   for node in nodeList:
      outstring += "%d %d 0 %f %f %f\n" % (node[0],node[0],node[1],node[2],node[3])
   outstring += "\nTri3 %d\n" % (len(triList))
   for tri in triList:
      outstring += "%d 0 %d %d %d\n" % (tri[0],tri[1],tri[2],tri[3])

   outfile = open(outfilename,"w")
   outfile.write(outstring)
   outfile.close()
   return

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]

   return nV

def reorderNodes(nodeList, triList, inPoint = [0.5,0.5,1.0]):
   for tid in range(len(triList)):
      id1 = triList[tid][1]
      id2 = triList[tid][2]
      id3 = triList[tid][3]
      p1 = [nodeList[id1][1], nodeList[id1][2], nodeList[id1][3]]
      p2 = [nodeList[id2][1], nodeList[id2][2], nodeList[id2][3]]
      p3 = [nodeList[id3][1], nodeList[id3][2], nodeList[id3][3]]

      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):
         triList[tid][1] = id3 
         triList[tid][3] = id1 

   return

def Medit2TriMesh (meshfile="corner.mesh", outfile="corner.lsm",inPoint = [0.5,0.5,1.0]):
   (nodeList, triList) = readMeditMesh(meshfile=meshfile)

   reorderNodes (nodeList, triList, inPoint = inPoint)

   writeTriMesh (nodeList,triList,outfilename=outfile)

   return

if __name__ == "__main__":
   if (len(argv)>=3):
      Medit2TriMesh(meshfile=argv[1],outfile=argv[2],inPoint=[0,10,0])

   else:
      Medit2TriMesh(meshfile="bucket.mesh",outfile="bucket.lsm",inPoint=[0,5,0])

GlossyBlue theme adapted by David Gilbert
Powered by PmWiki
www.000webhost.com