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

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

if __name__ == "__main__":
   pts,tris,bbox = readTriMeshFile(inputfile="Sample_0.msh")
   pts,tris,bbox = scale_and_translate(pts,tris,bbox,maxLength=10.0,centrePt=Vector3(0,0,0))
   writeTriMeshFile(pts,tris,"sample.msh")
GlossyBlue theme adapted by David Gilbert
Powered by PmWiki
www.000webhost.com