Re: [vmtk-users] branchsplit/summer school

Luca Antiga Mon, 09 May 2011 15:24:43 -0700

Dear Regine,

the results you are getting are very odd. This is not what it's supposed to be, and I can't reproduce it on tens of bifurcations, so it might be that there's something fishy with your particular model. Can you please send me the model so I can try it out first hand?

As for importing models from Gambit back to vmtk, there is no msh reader in vmtk (although there is a msh writer), but there is a GAMBIT neutral format reader in VTK (which is at the basis of vmtk).

If you want to try it out, just save your mesh in Gambit neutral format, and replace the vmtkmeshreader.py file in the lib directory in your vmtk installation with the one attached.

Use vmtkmeshreader -ifile foo.neu -f gambit --pipe vmtkmeshviewer -ofile foo.vtu

and see what comes out of it!

As for the Summer School, if you are interested in block structured hexahedral meshes we could arrange a session on it, for example based on vmtk/pyFormex.

Best regards

Luca

#!/usr/bin/env python

## Program:   VMTK
## Module:    $RCSfile: vmtkmeshreader.py,v $
## Language:  Python
## Date:      $Date: 2006/07/17 09:53:14 $
## Version:   $Revision: 1.16 $


##   Copyright (c) Luca Antiga, David Steinman. All rights reserved.
##   See LICENCE file for details.

##      This software is distributed WITHOUT ANY WARRANTY; without even 
##      the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR 
##      PURPOSE.  See the above copyright notices for more information.

import sys
import vtk
import vtkvmtk

import pypes

vmtkmeshreader = 'vmtkMeshReader'

class vmtkMeshReader(pypes.pypeScript):

    def __init__(self):

        pypes.pypeScript.__init__(self)

        self.Format = ''
        self.GuessFormat = 1
        self.InputFileName = ''
        self.Mesh = 0
        self.Output = 0
        
        self.GhostNodes = 1
        self.VolumeElementsOnly = 0

        self.CellEntityIdsArrayName = 'CellEntityIds'

        self.SetScriptName('vmtkmeshreader')
        self.SetScriptDoc('read a mesh and stores it in a vtkUnstructuredGrid object')
        self.SetInputMembers([
            ['Format','f','str',1,'["vtkxml","vtk","fdneut","ngneut","tecplot","tetgen","gambit"]','file format (fdneut - FIDAP neutral format, ngneut - Netgen neutral format)'],
            ['GuessFormat','guessformat','bool',1,'','guess file format from extension'],
            ['Mesh','i','vtkUnstructuredGrid',1,'','the input mesh'],
            ['InputFileName','ifile','str',1,'','input file name'],
            ['GhostNodes','ghostnodes','bool',1,'','store all nodes for 9-noded quads, 7-noded triangles, 27-noded hexahedra, 18-noded wedges; otherwise, store them as 8-noded quads, 6-noded triangles, 20-noded hexahedra, 15-noded wedges - fdneut only'],
            ['VolumeElementsOnly','volumeelementsonly','bool',1,'','only read volume elements - fdneut and ngneut'],
            ['CellEntityIdsArrayName','entityidsarray','str',1,'','name of the array where entity ids have to be stored - ngneut and tetgen']
            ])
        self.SetOutputMembers([
            ['Mesh','o','vtkUnstructuredGrid',1,'','the output mesh','vmtkmeshwriter'],
            ['CellEntityIdsArrayName','entityidsarray','str',1,'','name of the array where entity ids have been stored - ngneut and tetgen']
            ])

    def ReadTetGenMeshFile(self):
        if (self.InputFileName == ''):
            self.PrintError('Error: no InputFileName.')
        self.PrintLog('Reading TetGen mesh file.')
        import os.path
        inputFileName = os.path.splitext(self.InputFileName)[0]
        reader = vtkvmtk.vtkvmtkTetGenReader()
        reader.SetFileName(inputFileName)
        reader.SetBoundaryDataArrayName(self.CellEntityIdsArrayName)
        reader.Update()
        self.Mesh = reader.GetOutput()

    def ReadGAMBITMeshFile(self):
        if (self.InputFileName == ''):
            self.PrintError('Error: no InputFileName.')
        self.PrintLog('Reading GAMBIT mesh file.')
        reader = vtk.vtkGAMBITReader()
        reader.SetFileName(self.InputFileName)
        reader.Update()
        self.Mesh = reader.GetOutput()

    def ReadVTKMeshFile(self):
        if (self.InputFileName == ''):
            self.PrintError('Error: no InputFileName.')
        self.PrintLog('Reading VTK mesh file.')
        reader = vtk.vtkUnstructuredGridReader()
        reader.SetFileName(self.InputFileName)
        reader.Update()
        self.Mesh = reader.GetOutput()

    def ReadVTKXMLMeshFile(self):
        if (self.InputFileName == ''):
            self.PrintError('Error: no InputFileName.')
        self.PrintLog('Reading VTK XML mesh file.')
        reader = vtk.vtkXMLUnstructuredGridReader()
        reader.SetFileName(self.InputFileName)
        reader.Update()
        self.Mesh = reader.GetOutput()

    def ReadXdaMeshFile(self):
        if (self.InputFileName == ''):
            self.PrintError('Error: no InputFileName.')
        self.PrintError('Error: Xda reader not yet implemented.')
        return
#        self.PrintLog('Reading Xda mesh file.')
#        reader = vtkvmtk.vtkvmtkXdaReader()
#        reader.SetFileName(self.InputFileName)
#        reader.Update()
#        self.Mesh = reader.GetOutput()

    def ReadFDNEUTMeshFile(self):
        if (self.InputFileName == ''):
            self.PrintError('Error: no InputFileName.')
        self.PrintLog('Reading FDNEUT mesh file.')
        reader = vtkvmtk.vtkvmtkFDNEUTReader()
        reader.SetFileName(self.InputFileName)
        reader.SetGhostNodes(self.GhostNodes)
        reader.SetSingleCellDataEntityArrayName(self.CellEntityIdsArrayName)
        reader.SetVolumeElementsOnly(self.VolumeElementsOnly)
        reader.Update()
        self.Mesh = reader.GetOutput()

    def ReadNGNEUTMeshFile(self):
        if (self.InputFileName == ''):
            self.PrintError('Error: no InputFileName.')
        self.PrintLog('Reading ngneut mesh file.')
        f=open(self.InputFileName, 'r')
        self.Mesh = vtk.vtkUnstructuredGrid()
        self.MeshPoints = vtk.vtkPoints()
        line = f.readline()
        numberOfPoints = int(line)
        self.MeshPoints.SetNumberOfPoints(numberOfPoints)
        for i in range(numberOfPoints):
            line = f.readline()
            splitLine = line.strip().split()
            point = [float(splitLine[0]),float(splitLine[1]),float(splitLine[2])]
            self.MeshPoints.SetPoint(i,point)
        self.Mesh.SetPoints(self.MeshPoints)
        self.Mesh.Allocate(numberOfPoints*4,1000)
        line = f.readline()
        numberOfVolumeCells = int(line)
        cellEntityIdArray = vtk.vtkIntArray()
        cellEntityIdArray.SetName(self.CellEntityIdsArrayName)
        for i in range(numberOfVolumeCells):
            line = f.readline()
            splitLine = line.strip().split()
            cellType = -1
            numberOfCellPoints = len(splitLine)-1
            pointIds = vtk.vtkIdList()
            if numberOfCellPoints == 4:
                cellType = 10
#                pointIds.InsertNextId(int(splitLine[0+1])-1)
#                pointIds.InsertNextId(int(splitLine[1+1])-1)
#                pointIds.InsertNextId(int(splitLine[2+1])-1)
#                pointIds.InsertNextId(int(splitLine[3+1])-1)
                pointIds.InsertNextId(int(splitLine[1+1])-1)
                pointIds.InsertNextId(int(splitLine[2+1])-1)
                pointIds.InsertNextId(int(splitLine[3+1])-1)
                pointIds.InsertNextId(int(splitLine[0+1])-1)
            elif numberOfCellPoints == 10:
                cellType = 24
#                pointIds.InsertNextId(int(splitLine[0+1])-1)
#                pointIds.InsertNextId(int(splitLine[1+1])-1)
#                pointIds.InsertNextId(int(splitLine[2+1])-1)
#                pointIds.InsertNextId(int(splitLine[3+1])-1)
#                pointIds.InsertNextId(int(splitLine[4+1])-1)
#                pointIds.InsertNextId(int(splitLine[7+1])-1)
#                pointIds.InsertNextId(int(splitLine[5+1])-1)
#                pointIds.InsertNextId(int(splitLine[6+1])-1)
#                pointIds.InsertNextId(int(splitLine[8+1])-1)
#                pointIds.InsertNextId(int(splitLine[9+1])-1)
                pointIds.InsertNextId(int(splitLine[1+1])-1)
                pointIds.InsertNextId(int(splitLine[2+1])-1)
                pointIds.InsertNextId(int(splitLine[3+1])-1)
                pointIds.InsertNextId(int(splitLine[0+1])-1)
                pointIds.InsertNextId(int(splitLine[7+1])-1)
                pointIds.InsertNextId(int(splitLine[9+1])-1)
                pointIds.InsertNextId(int(splitLine[8+1])-1)
                pointIds.InsertNextId(int(splitLine[4+1])-1)
                pointIds.InsertNextId(int(splitLine[5+1])-1)
                pointIds.InsertNextId(int(splitLine[6+1])-1)
            entityId = int(splitLine[0])
            cellEntityIdArray.InsertNextValue(entityId)
            self.Mesh.InsertNextCell(cellType,pointIds)
        if self.VolumeElementsOnly == 0:
            line = f.readline()
            numberOfSurfaceCells = int(line)
            for i in range(numberOfSurfaceCells):
                line = f.readline()
                splitLine = line.strip().split()
                cellType = -1
                numberOfCellPoints = len(splitLine)-1
                pointIds = vtk.vtkIdList()
                if numberOfCellPoints == 3:
                    cellType = 5
                    for j in range(numberOfCellPoints):
                        pointIds.InsertNextId(int(splitLine[j+1])-1)
                elif numberOfCellPoints == 6:
                    cellType = 22
                    for j in range(3):
                        pointIds.InsertNextId(int(splitLine[j+1])-1)
                    pointIds.InsertNextId(int(splitLine[5+1])-1)
                    pointIds.InsertNextId(int(splitLine[3+1])-1)
                    pointIds.InsertNextId(int(splitLine[4+1])-1)
                entityId = int(splitLine[0])
                cellEntityIdArray.InsertNextValue(entityId)
                self.Mesh.InsertNextCell(cellType,pointIds)
        self.Mesh.Squeeze()
        self.Mesh.GetCellData().AddArray(cellEntityIdArray)

    def ReadTecplotMeshFile(self):
        self.PrintLog('Reading Tecplot surface file.')
        if self.InputFileName[-3:] == '.gz':
            import gzip
            f=gzip.open(self.InputFileName, 'r')
        else:
            f=open(self.InputFileName, 'r')
        line = f.readline()
        if line.split()[0] == 'TITLE':
            line = f.readline()
        if (line.split()[0] == 'VARIABLES') | (line.split('=')[0] == 'VARIABLES'):
            arrayNames = line.split('=')[1].strip().split(',')
            arrayNames[0:3] = []
            self.PrintLog("ArrayNames" + str(arrayNames))
            line = f.readline()
        if line.split()[0] == 'ZONE':
            lineNid = line.find('N=')
            lineN = line[lineNid : lineNid+line[lineNid:].find(',') ].split('=')[1]
            numberOfNodes = int(lineN)
            lineEid = line.find('E=')
            lineE = line[lineEid : lineEid+line[lineEid:].find(',') ].split('=')[1]
            numberOfElements = int(lineE)
        self.PrintLog("Reading " + str(numberOfNodes)+" nodes.")
        points = vtk.vtkPoints()
        cells = vtk.vtkCellArray()
        points.SetNumberOfPoints(numberOfNodes)
        self.Mesh = vtk.vtkUnstructuredGrid()
        self.Mesh.SetPoints(points)
        for arrayName in arrayNames:
            array = vtk.vtkDoubleArray()
            array.SetName(arrayName)
            array.SetNumberOfTuples(numberOfNodes)
            self.Mesh.GetPointData().AddArray(array)
        data = f.read().split()
        dataCounter = 0
        for i in range(numberOfNodes):
            point = [float(data[dataCounter]),float(data[dataCounter+1]),float(data[dataCounter+2])]
            dataCounter += 3
            points.SetPoint(i,point)
            for j in range(len(arrayNames)):
                self.Mesh.GetPointData().GetArray(arrayNames[j]).SetComponent(i,0,float(data[dataCounter]))
                dataCounter += 1
        self.PrintLog("Reading " + str(numberOfElements)+" elements.")
        cellIds = vtk.vtkIdList()
        for i in range(numberOfElements):
            cellIds.Initialize()
            cellIds.InsertNextId(int(data[dataCounter])-1)
            dataCounter += 1
            cellIds.InsertNextId(int(data[dataCounter])-1)
            dataCounter += 1
            cellIds.InsertNextId(int(data[dataCounter])-1)
            dataCounter += 1
            cellIds.InsertNextId(int(data[dataCounter])-1)
            dataCounter += 1
            legal = 1
            for j in range(cellIds.GetNumberOfIds()):
                if cellIds.GetId(j) < 0:
                    legal = 0
                    break
            if not legal:
                continue
            cells.InsertNextCell(cellIds)
        self.Mesh.SetCells(10,cells)
        self.Mesh.Update()

    def Execute(self):

        extensionFormats = {'vtu':'vtkxml', 
                            'vtkxml':'vtkxml', 
                            'vtk':'vtk',
                            'FDNEUT':'fdneut',
                            'xda':'xda',
                            'neu':'ngneut',
                            'gneu':'gambit',
                            'tec':'tecplot',
                            'node':'tetgen',
                            'ele':'tetgen'}

        if self.InputFileName == 'BROWSER':
            import tkFileDialog
            import os.path
            initialDir = pypes.pypeScript.lastVisitedPath
            self.InputFileName = tkFileDialog.askopenfilename(title="Input mesh",initialdir=initialDir)
            pypes.pypeScript.lastVisitedPath = os.path.dirname(self.InputFileName)
            if not self.InputFileName:
                self.PrintError('Error: no InputFileName.')

        if self.GuessFormat and self.InputFileName and not self.Format:
            import os.path
            extension = os.path.splitext(self.InputFileName)[1]
            if extension:
                extension = extension[1:]
                if extension in extensionFormats.keys():
                    self.Format = extensionFormats[extension]

        if (self.Format == 'vtk'):
            self.ReadVTKMeshFile()
        elif (self.Format == 'vtkxml'):
            self.ReadVTKXMLMeshFile()
        elif (self.Format == 'gambit'):
            self.ReadGAMBITMeshFile()
        elif (self.Format == 'fdneut'):
            self.ReadFDNEUTMeshFile()
        elif (self.Format == 'ngneut'):
            self.ReadNGNEUTMeshFile()
        elif (self.Format == 'xda'):
            self.ReadXdaMeshFile()
        elif (self.Format == 'tecplot'):
            self.ReadTecplotMeshFile()
        elif (self.Format == 'tetgen'):
            self.ReadTetGenMeshFile()
        else:
            self.PrintError('Error: unsupported format '+ self.Format + '.')

        if self.Mesh.GetSource():
            self.Mesh.GetSource().UnRegisterAllOutputs()

        self.Output = self.Mesh


if __name__=='__main__':
    main = pypes.pypeMain()
    main.Arguments = sys.argv
    main.Execute()

On May 3, 2011, at 12:48 PM, richschm...@web.de wrote:

Dear Luca,

thank you for your help!

Unfortunately I always get strange results, if I place the source seed at the inlet (I added some pics: Branchsplit1,2,3). The results of the branch splitting vary as well, if I change the order of the seed point positioning (pic: Branchsplit4). Therefore I started to place them at the bifurcation apex. Is there any special trick?

You said I could merge the branches afterwards again. I meshed my branches using Gambit and saved the result as a *.msh file (unfortunately there are not that many possible data types I can use, because part of my geometry is “non-real” and therefore it vanishes, if I save it as *.step or *.iges file). Is there any possibility to import such a *.msh file at vmtk? I tried to do this and failed several times and it’s also pretty hard to find a freeware to convert the *.msh file into *.stl again. Is it possible to merge these meshed branches afterwards anyway?

The topics of the summer school sound quite interesting. Will one of the lectures or tutorials also include the generation of block structured hexahedral grids in a vessel geometry?

Thanks a million in advance!
With best regards,
Regine
Von: "Luca Antiga" <luca.ant...@gmail.com>
Gesendet: 07.04.2011 11:55:10
An: richschm...@web.de
Betreff: Re: [vmtk-users] blockstructured hexahedral mesh at bifurcation

Dear Regine,
The pipe is correct the way you write it.
The problem with branch splitting the way you are doing it is that you should place your source seed down at the CCA inlet, not at the bifurcation apex.
Don't worry if the centerlines do not meet at the bifurcation "center" (whatever that is), it's correct that way (you can always merge them *after splitting the surface* using vmtkcenterlinemerge).

If you want the cut surfaces to match exactly, you have to increase the density of your surface using vmtksurfacesubdivision, using 1 or 2 levels of subdivision. This way the cut will be more "precise" (up to a numeric tolerance).

Let me know if this answers your question, best regards

Luca

On Apr 6, 2011, at 1:50 PM, richschm...@web.de wrote:

Dear Luca,

thank you for your help! Sounds good. I will try my best with pyformex!

I have two more questions dealing with vmtk branch splitting. I hope you (or somebody else) can help me.

When I was splitting all three branches of my bifurcation I never managed to start from the same centerline generation. I always had to fix the red dots for the source points newly and therefore the cut surfaces differ a bit, that means the branches don’t really match at the bifurcation. I was using this code:

“vmtksurfacereader -ifile E:\Regine\ Model_clipped.stl --pipe vmtkcenterlines --pipe vmtkbranchextractor --pipe vmtkbranchclipper -groupids 2 -ofile E:\Regine\Model_branchextract1.stl
vmtksurfaceviewer -ifile E:\Regine \Model_branchextract1.stl -array GroupIds”

I tried different variations (for the single group ids) to use the same starting point for the branch splitting, but failed.

The other thing is concerning the position of my branch splitting. It doesn’t look like the one shown at vmtk.org and I am concerned that my mesh won't be good therefore. I tried different positions of the seed point, but it didn't work. I attached two pics showing where I was positioning the seed points and how my branchsplit looks afterwards.

It would be really nice, if you could help!

With best regards,
Regine

Von: "Luca Antiga" <luca.ant...@gmail.com>
Gesendet: 05.04.2011 08:53:16
An: richschm...@web.de
Betreff: Re: [vmtk-users] blockstructured hexahedral mesh at bifurcation

Dear Regine,
if you are interested in a block-structured hexahedral mesh, I strongly suggest you check out pyformex, which does exactly what you look for, http://pyformex.berlios.de/.
You can import an stl file and generate an all-hex mesh.
Hope this helps.

Luca

On Apr 3, 2011, at 10:50 PM, richschm...@web.de wrote:

Dear all,

I am using the branch splitting to seperate the branches of a vessel bifurcation. My plan is to create a blockstructured hexahedral mesh for such a bifurcation. Since this appeared to be not that easy considering the whole bifurcation (I failed several times), I hoped this would maybe get easier, if I split the branches. Unfortunately I am not that experienced with mesh generation and cfd in general. Does any of you have experience how I can proceed now to get such a blockstructured hexahedral mesh? The only commercial software I have access to are Gambit and TGrid. It would be really nice, if you could help me!

Thanks a million in advance!
Regine

Empfehlen Sie WEB.DE DSL Ihren Freunden und Bekannten und wir
belohnen Sie mit bis zu 50,- Euro! https://freundschaftswerbung.web.de
------------------------------------------------------------------------------
Create and publish websites with WebMatrix
Use the most popular FREE web apps or write code yourself;
WebMatrix provides all the features you need to develop and
publish your website. http://p.sf.net/sfu/ms-webmatrix-sf
_______________________________________________
vmtk-users mailing list
vmtk-users@lists.sourceforge.net
https://lists.sourceforge.net/lists/listinfo/vmtk-users


WEB.DE DSL Doppel-Flat ab 19,99 €/mtl.! Jetzt mit
gratis Handy-Flat! http://produkte.web.de/go/DSL_Doppel_Flatrate/2
<SeedPoints_Centerline.jpg><Branchsplit.jpg>


Schon gehört? WEB.DE hat einen genialen Phishing-Filter in die
Toolbar eingebaut! http://produkte.web.de/go/toolbar

<Branchsplit1.jpg><Branchsplit2.jpg><Branchsplit3.jpg><Branchsplit4.jpg>

------------------------------------------------------------------------------
Achieve unprecedented app performance and reliability
What every C/C++ and Fortran developer should know.
Learn how Intel has extended the reach of its next-generation tools
to help boost performance applications - inlcuding clusters.
http://p.sf.net/sfu/intel-dev2devmay

_______________________________________________
vmtk-users mailing list
vmtk-users@lists.sourceforge.net
https://lists.sourceforge.net/lists/listinfo/vmtk-users

Re: [vmtk-users] branchsplit/summer school

Reply via email to