convert2TetGen.py

Go to the documentation of this file.
00001 # (c) 2010 LGPL
00002 
00003 #Make mesh of pn junction in TetGen format
00004 import FreeCAD, FreeCADGui, Part, Mesh
00005 App = FreeCAD # shortcut
00006 Gui = FreeCADGui # shortcut
00007 
00008 def exportMeshToTetGenPoly(meshToExport,filePath,beVerbose=1):
00009     """Export mesh to  TetGen *.poly file format"""
00010     ## Part 1 - write node list to output file
00011     if beVerbose == 1:
00012             FreeCAD.Console.PrintMessage("\nExport of mesh to TetGen file ...")
00013     (allVertices,allFacets) = meshToExport.Topology
00014     f = open(filePath, 'w')
00015     f.write("# This file was generated from FreeCAD geometry\n")
00016     f.write("# Part 1 - node list\n")
00017     f.write("%(TotalNumOfPoints)i  %(NumOfDimensions)i  %(NumOfProperties)i  %(BoundaryMarkerExists)i\n" % \
00018             {'TotalNumOfPoints':len(allVertices), \
00019              'NumOfDimensions':3, \
00020              'NumOfProperties':0, \
00021              'BoundaryMarkerExists':0})
00022     for PointIndex in range(len(allVertices)):
00023         f.write("%(PointIndex)5i %(x) e %(y) e %(z) e\n" % \
00024                 {'PointIndex':PointIndex, \
00025                 'x':allVertices[PointIndex].x, \
00026                 'y':allVertices[PointIndex].y, \
00027                 'z':allVertices[PointIndex].z})
00028 
00029     ## Find out BoundaryMarker for each facet. If edge connects only two facets,
00030     # then this facets should have the same BoundaryMarker
00031     BoundaryMarkerExists = 1
00032     PointList=[allFacets[0][1],allFacets[0][0]]
00033     PointList.sort()
00034     EdgeFacets = {(PointList[0],PointList[1]):set([0])}
00035     Edge = []
00036 
00037     # Finde all facets for each edge
00038     for FacetIndex in range(len(allFacets)):
00039         Facet = allFacets[FacetIndex]
00040         for i in range(0,-len(Facet),-1):
00041             tmpEdge = [Facet[i],Facet[i+1]]
00042             tmpEdge.sort()
00043             Edge.append(tmpEdge)
00044         for i in range(len(Edge)):
00045             EdgeIndex = (Edge[i][0],Edge[i][1])
00046             if EdgeIndex in EdgeFacets:
00047                 EdgeFacets[EdgeIndex].add(FacetIndex)               
00048             else:
00049                 EdgeFacets[EdgeIndex] = set([FacetIndex])                
00050         Edge = []
00051     
00052     # Find BoundaryMarker for each facet
00053     BoundaryMarker = []
00054     for index in range(len(allFacets)):
00055         BoundaryMarker.append(0)
00056     MinMarker = -1
00057     InitialFacet = 0
00058     BoundaryMarker[InitialFacet] = MinMarker
00059     EdgeKeys = EdgeFacets.keys()
00060     disconnectedEdges = len(EdgeKeys)
00061     if beVerbose == 1:
00062         FreeCAD.Console.PrintMessage('\nBoundaryMarker:'+repr(BoundaryMarker)+' '+repr(len(EdgeFacets)))
00063     searchForPair = 1
00064 
00065     # Main loop: first search for all complementary facets, then fill one branch and repeat while edges are available
00066     while len(EdgeFacets)>0:        
00067         removeEdge = 0
00068         for EdgeIndex in EdgeKeys:
00069             if len(EdgeFacets[EdgeIndex]) == 1:
00070                 removeEdge = 1
00071                 break
00072             if len(EdgeFacets[EdgeIndex]) == 2:
00073                 FacetPair = []
00074                 for facet in EdgeFacets[EdgeIndex]:
00075                     FacetPair.append(facet)
00076                 if (BoundaryMarker[FacetPair[0]]==0) and (BoundaryMarker[FacetPair[1]]==0):
00077                     continue
00078                 if (BoundaryMarker[FacetPair[0]]!=0) and (BoundaryMarker[FacetPair[1]]!=0):
00079                     removeEdge = 1
00080                     break
00081                 if (BoundaryMarker[FacetPair[0]]!=0):
00082                     BoundaryMarker[FacetPair[1]] = BoundaryMarker[FacetPair[0]]
00083                 else:
00084                     BoundaryMarker[FacetPair[0]] = BoundaryMarker[FacetPair[1]]
00085                 removeEdge = 1
00086                 break
00087             if searchForPair == 1:
00088                 continue
00089             FacetTree = []
00090             AllMarkers = 1
00091             MarkerSum = 0
00092             for facet in EdgeFacets[EdgeIndex]:
00093                 FacetTree.append(facet)
00094                 MarkerSum += BoundaryMarker[facet]
00095             if MarkerSum == 0:
00096                 continue
00097             for facet in EdgeFacets[EdgeIndex]:
00098                 if BoundaryMarker[facet] == 0:
00099                     MinMarker -= 1
00100                     BoundaryMarker[facet] = MinMarker
00101                     searchForPair = 1
00102             removeEdge = 1
00103             break
00104         if removeEdge == 1:
00105             del EdgeFacets[EdgeIndex]
00106             EdgeKeys = EdgeFacets.keys()
00107             continue
00108         searchForPair = 0
00109     # End of main loop
00110     if beVerbose == 1:
00111         FreeCAD.Console.PrintMessage('\nNew BoundaryMarker:'+repr(BoundaryMarker)+' '+repr(len(EdgeFacets)))
00112 
00113     ## Part 2 - write all facets to *.poly file
00114     f.write("# Part 2 - facet list\n")
00115     f.write("%(TotalNumOfFacets)i  %(BoundaryMarkerExists)i\n" %\
00116             {'TotalNumOfFacets':len(allFacets),\
00117              'BoundaryMarkerExists':BoundaryMarkerExists})
00118     for FacetIndex in range(len(allFacets)):
00119         f.write("# FacetIndex = %(Index)i\n"%{'Index':FacetIndex})
00120         f.write("%(NumOfPolygons)3i "%{'NumOfPolygons':1})
00121         if BoundaryMarkerExists == 1:
00122             f.write("0 %(BoundaryMarker)i"%{'BoundaryMarker':BoundaryMarker[FacetIndex]})
00123         f.write("\n%(NumOfConers)3i  "%{'NumOfConers':len(allFacets[FacetIndex])})
00124         for PointIndex in range(len(allFacets[FacetIndex])):
00125             #        f.write(repr(allFacets[FacetIndex][PointIndex]))
00126             f.write("%(PointIndex)i "%{'PointIndex':allFacets[FacetIndex][PointIndex]})
00127         f.write("\n")
00128     ## Part 3 and Part 4 are zero
00129     f.write("# Part 3 - the hole list.\n# There is no hole in bar.\n0\n")
00130     f.write("# Part 4 - the region list.\n# There is no region defined.\n0\n")
00131     f.write("# This file was generated from FreeCAD geometry\n")
00132     f.close()
00133 
00134 
00135 def export(objectslist,filename):
00136     """Called when freecad exports a mesh to poly format"""
00137     for obj in objectslist:
00138         if isinstance(obj,Mesh.Feature):
00139             exportMeshToTetGenPoly(obj.Mesh,filename,False)
00140             break
00141 
00142 def createMesh():
00143     ## ========================  Script beginning...  ========================
00144     beVerbose = 1
00145     if beVerbose == 1:
00146         FreeCAD.Console.PrintMessage("\n\n\n\n\n\n\n\nScript starts...") 
00147     ## Geometry defenition
00148     # Define objects names
00149     PyDocumentName = "pnJunction"
00150     PSideBoxName = "PSide"
00151     NSideBoxName = "NSide"
00152     DepletionBoxName = "Depletion"
00153     SurfDepletionBoxName = "SurfDepletion"
00154     OxideBoxName = "Oxide"
00155     AdsorbtionBoxName = "Adsorbtion"
00156     pnMeshName = "pnMesh"
00157 
00158     # Init objects
00159     if beVerbose == 1:
00160         FreeCAD.Console.PrintMessage("\nInit Objects...") 
00161     # App.closeDocument(App.ActiveDocument.Label) #closeDocument after restart of macro. Needs any ActiveDocument.
00162     AppPyDoc = App.newDocument(PyDocumentName)
00163     NSideBox = AppPyDoc.addObject("Part::Box",NSideBoxName)
00164     PSideBox = AppPyDoc.addObject("Part::Box",PSideBoxName)
00165     DepletionBox = AppPyDoc.addObject("Part::Box",DepletionBoxName)
00166     SurfDepletionBox = AppPyDoc.addObject("Part::Box",SurfDepletionBoxName)
00167     OxideBox = AppPyDoc.addObject("Part::Box",OxideBoxName)
00168     AdsorbtionBox = AppPyDoc.addObject("Part::Box",AdsorbtionBoxName)
00169     pnMesh = AppPyDoc.addObject("Mesh::Feature",pnMeshName)
00170 
00171     BoxList=[NSideBox, DepletionBox, PSideBox,  OxideBox, AdsorbtionBox, SurfDepletionBox]
00172     NSideBoxMesh = Mesh.Mesh()
00173     PSideBoxMesh = Mesh.Mesh()
00174     DepletionBoxMesh = Mesh.Mesh()
00175     SurfDepletionBoxMesh = Mesh.Mesh()
00176     OxideBoxMesh = Mesh.Mesh()
00177     AdsorbtionBoxMesh = Mesh.Mesh()
00178     BoxMeshList = [NSideBoxMesh, DepletionBoxMesh, PSideBoxMesh, OxideBoxMesh, AdsorbtionBoxMesh, SurfDepletionBoxMesh]
00179     if beVerbose == 1:
00180         if len(BoxList)!=len(BoxMeshList):
00181             FreeCAD.Console.PrintMessage("\n ERROR! Input len() of BoxList and BoxMeshList is not the same! ")
00182 
00183     ## Set sizes in nanometers
00184     if beVerbose == 1:
00185         FreeCAD.Console.PrintMessage("\nSet sizes...") 
00186     tessellationTollerance = 0.05
00187     ModelWidth = 300
00188     BulkHeight = 300
00189     BulkLength = 300
00190     DepletionSize = 50
00191     OxideThickness = 5
00192     AdsorbtionThickness = 10
00193 
00194     # Big volumes of n and p material
00195     NSideBox.Height = BulkHeight #Z-direction
00196     NSideBox.Width = ModelWidth  #Y-direction = const
00197     NSideBox.Length = BulkLength  #X-direction
00198     PSideBox.Height = BulkHeight
00199     PSideBox.Width = ModelWidth
00200     PSideBox.Length = BulkLength
00201     # Thin depletion layer between
00202     DepletionBox.Height = BulkHeight
00203     DepletionBox.Width = ModelWidth
00204     DepletionBox.Length = DepletionSize*2
00205     # Surface deplation layer
00206     SurfDepletionBox.Height = DepletionSize
00207     SurfDepletionBox.Width = ModelWidth
00208     SurfDepletionBox.Length = BulkLength*2 + DepletionSize*2
00209     # Oxide on the top
00210     OxideBox.Height = OxideThickness
00211     OxideBox.Width = ModelWidth
00212     OxideBox.Length = BulkLength*2 + DepletionSize*2
00213     # Adsorbtion layer
00214     AdsorbtionBox.Height = AdsorbtionThickness
00215     AdsorbtionBox.Width = ModelWidth
00216     AdsorbtionBox.Length = BulkLength*2 + DepletionSize*2
00217 
00218     # Object placement
00219     Rot = App.Rotation(0,0,0,1)
00220     NSideBox.Placement = App.Placement(App.Vector(0,0,-BulkHeight),Rot)
00221     PSideBox.Placement = App.Placement(App.Vector(DepletionSize*2+BulkLength,0,-BulkHeight),Rot)
00222     DepletionBox.Placement = App.Placement(App.Vector(BulkLength,0,-BulkHeight),Rot)
00223     SurfDepletionBox.Placement = App.Placement(App.Vector(0,0,0),Rot)
00224     OxideBox.Placement = App.Placement(App.Vector(0,0,DepletionSize),Rot)
00225     AdsorbtionBox.Placement = App.Placement(App.Vector(0,0,DepletionSize+OxideThickness),Rot)
00226 
00227     ## Unite
00228     if beVerbose == 1:
00229         FreeCAD.Console.PrintMessage("\nFuse objects...") 
00230     fuseShape = BoxList[0].Shape
00231     for index in range(1,len(BoxList),1):
00232         fuseShape = fuseShape.fuse(BoxList[index].Shape)
00233     nmesh = Mesh.Mesh()
00234     nmesh.addFacets(fuseShape.tessellate(tessellationTollerance))
00235 
00236     # for index in range(len(BoxList)): 
00237     for index in range(len(BoxList)-1): # Manual hack
00238         BoxMeshList[index].addFacets(BoxList[index].Shape.tessellate(tessellationTollerance))
00239         nmesh.addMesh(BoxMeshList[index])
00240 
00241     nmesh.removeDuplicatedPoints()
00242     nmesh.removeDuplicatedFacets()
00243     pnMesh.Mesh = nmesh
00244 
00245     # Hide all boxes
00246     for box in BoxList:
00247         Gui.hideObject(box)
00248     # # Remove all boxes
00249     # for box in BoxList:
00250     #     App.ActiveDocument.removeObject(box.Name)
00251     
00252     # Update document
00253     AppPyDoc.recompute()
00254 
00255     ## export to TenGen *.poly (use File|Export instead)
00256     #filePath = "/home/tig/tmp/tetgen/pnJunction.poly"
00257     #exportMeshToTetGenPoly(pnMesh.Mesh,filePath,beVerbose)
00258 
00259     Gui.activeDocument().activeView().viewAxometric()
00260     Gui.SendMsgToActiveView("ViewFit")
00261     if beVerbose == 1:
00262         FreeCAD.Console.PrintMessage("\nScript finished without errors.")