[Insight-users] itk::fem
Luis Ibanez
luis.ibanez at kitware.com
Wed Nov 21 16:00:09 EST 2007
Hi Vincent,
Thanks for sharing these files.
A suggestion:
It will be great to post them as a submission
to the Insight Journal. The will be useful to
many users.
Regards,
Luis
-------------------------
Vincent A. Magnotta wrote:
> Markus,
>
> We have implemented something similar using the itk::fem framework. We
> have always treated the deforming mesh as a solid meaning that it
> contains hex or tet elements that compose the structure. In this case we
> will load a vtkUnstructured grid convert it to an ITK mesh and then to a
> itk::fem mesh. We then deform the mesh by applying loads to the surface
> nodes of the mesh while applying a boundary condition to the center node
> to hold it fixed.
>
> I have attached code for the conversion from vtk into an itk::fem mesh.
>
>
> Vince
>
>
>
> On Tue, 2007-11-20 at 10:12 +0100, Markus Mehrwald wrote:
>
>>Sorry...the whole Mail again!
>>
>>Hello everyone,
>>
>>I am looking for an example for the itk:fem package. Following a short
>>description of of my problem.
>>I have a sphere as VTK Polydata. This sphere should be transfered into a
>>image which I can use with itk:fem (actually I do not understand why I
>>cannot use a mesh). To the center of the sphere should a force be
>>applied in the outward direction that the sphere is growing. The bigger
>>sphere should be transformed back into a VTK Polydata.
>>Has anyone an example for this or even for parts of the whole problem.
>>
>>Thank you very much in advance,
>>Markus
>>
>>
>>
>>------------------------------------------------------------------------
>>
>>#ifndef __ItkMeshToFEMMesh_h__
>>#define __ItkMeshToFEMMesh_h__
>>
>>#include "itkLightProcessObject.h"
>>#include <iostream>
>>#include "itkMesh.h"
>>#include "itkAutomaticTopologyMeshSource.h"
>>#include "itkTriangleCell.h"
>>#include "itkQuadrilateralCell.h"
>>#include "itkTetrahedronCell.h"
>>#include "itkHexahedronCell.h"
>>#include "itkFEM.h"
>>#include "itkCellInterface.h"
>>
>>namespace itk
>>{
>>
>>template<class TInputMesh>
>>class ITK_EXPORT ItkMeshToFEMMesh : public LightProcessObject
>> {
>>
>> public:
>>
>> typedef ItkMeshToFEMMesh Self;
>> typedef SmartPointer<Self> Pointer;
>>
>> typedef TInputMesh MeshType;
>> typedef typename MeshType::Pointer MeshTypePointer;
>> typedef typename MeshType::CellType CellType;
>> typedef typename MeshType::PointType PointType;
>> typedef typename MeshType::CellPixelType CellPixelType;
>> typedef typename MeshType::CellsContainerPointer CellsContainerPointer;
>> typedef typename MeshType::CellsContainer::ConstIterator CellIterator;
>> typedef typename itk::VertexCell< CellType > VertexType;
>> typedef typename itk::HexahedronCell< CellType > HexahedronType;
>> typedef typename HexahedronType::PointIdIterator PointIdIterator;
>>
>> typedef itk::fem::Solver SolverType;
>> typedef SolverType* SolverPointerType;
>> typedef itk::fem::MaterialLinearElasticity MaterialType;
>> typedef itk::fem::Element3DC0LinearHexahedronStrain HexElementType;
>> typedef HexElementType::Node HexNodeType;
>> typedef itk::fem::Element2DC0LinearTriangularStrain TriangElementType;
>> typedef TriangElementType::Node TriangNodeType;
>>
>> /**
>> Useful Internal Typedefs
>>
>> typedef itk::AutomaticTopologyMeshSource< MeshType > MeshSourceType;
>> typedef MeshSourceType::Pointer MeshSourceTypePointer;
>>*/
>> itkNewMacro( Self );
>>
>> itkSetObjectMacro ( Input, MeshType );
>>
>> itkSetStringMacro( FileName );
>> void Update ( );
>> SolverPointerType GetOutput ( );
>>
>>protected:
>>
>> MeshTypePointer m_Input;
>> std::string m_FileName;
>> SolverPointerType m_Solver ;
>>
>> ItkMeshToFEMMesh( );
>> virtual ~ItkMeshToFEMMesh( );
>> };
>>}
>>
>>#ifndef ITK_MANUAL_INSTANTIATION
>>#include "ItkMeshToFEMMesh.txx"
>>#endif
>>
>>#endif
>>
>>
>>------------------------------------------------------------------------
>>
>>#ifndef __ItkMeshToFEMMesh_cxx
>>#define __ItkMeshToFEMMesh_cxx
>>
>>#include <iostream>
>>#include "ItkMeshToFEMMesh.h"
>>
>>namespace itk
>>{
>>
>>template <class TInputMesh>
>>ItkMeshToFEMMesh<TInputMesh>
>>::ItkMeshToFEMMesh()
>> {
>> m_FileName = "";
>> m_Solver = new SolverType;
>> }
>>
>>template <class TInputMesh>
>>ItkMeshToFEMMesh<TInputMesh>
>>::~ItkMeshToFEMMesh()
>> {
>> delete m_Solver;
>> }
>>
>>template <class TInputMesh>
>>typename ItkMeshToFEMMesh<TInputMesh>::SolverPointerType
>>ItkMeshToFEMMesh<TInputMesh>
>>::GetOutput( )
>> {
>> return m_Solver;
>> }
>>
>>
>>template <class TInputMesh>
>>void
>>ItkMeshToFEMMesh<TInputMesh>
>>::Update( )
>>{
>> // Create dummy material
>> MaterialType::Pointer mat = MaterialType::New();
>> mat->GN = 0;
>> mat->E = 200E6;
>> mat->A = 1.0;
>> mat->h = 1.0;
>> mat->I = 1.0;
>> mat->nu = 0.2;
>> mat->RhoC = 1.0;
>>
>> // Create element type
>> HexElementType::Pointer hexElement = HexElementType::New( );
>> HexElementType::Pointer hexElement2;
>> hexElement->m_mat = dynamic_cast<MaterialType*>( mat );
>>
>> TriangElementType::Pointer triangElement = TriangElementType::New( );
>> TriangElementType::Pointer triangElement2;
>> triangElement->m_mat = dynamic_cast<MaterialType*>( mat );
>>
>> // Convert mesh points into nodes
>> PointType* ptr;
>> PointType pt;
>> ptr = &pt;
>>
>> int numofpoints = m_Input->GetNumberOfPoints();
>> //std::cerr << "# points: " << numofpoints << std::endl;
>>
>>
>> // Convert cells into elements
>> CellsContainerPointer cellList = m_Input->GetCells();
>> CellIterator cells = cellList->Begin();
>>
>>
>> //std::cerr << m_Input << std::endl;
>> int numberOfCells = m_Input->GetNumberOfCells( );
>> //std::cerr << "# cells: " << numberOfCells << std::endl;
>>
>> bool Flag = true;
>>
>> for (int k=0; k < numberOfCells; k++)
>> {
>> CellType* cellPtr = cells.Value();
>> //std::cout << "CELL TYPE = " << cellPtr->GetType() << " Hex Element Type " <<
>> //CellType::HEXAHEDRON_CELL << std::endl;
>>
>> switch( cellPtr->GetType( ) )
>> {
>> case CellType::HEXAHEDRON_CELL:
>> {
>> //std::cout << "Add Hex Element" << std::endl;
>> if ( Flag ) // To make sure that the nodes are created just once
>> {
>> for (int j=0; j < numofpoints; j++)
>> {
>> m_Input->GetPoint(j, ptr);
>>
>> HexNodeType::Pointer hexNode;
>> hexNode = new HexNodeType(pt[0], pt[1], pt[2]);
>> hexNode->GN = j;
>> m_Solver->node.push_back(itk::fem::FEMP<HexNodeType>(hexNode));
>> }
>> Flag = false;
>> }
>> PointIdIterator pointIt = cellPtr->PointIdsBegin() ;
>> int i=0;
>> hexElement2 = dynamic_cast<HexElementType*> (hexElement->Clone());
>> while (pointIt != cellPtr->PointIdsEnd() )
>> {
>> hexElement2->SetNode(i, m_Solver->node.Find( *pointIt ));
>> pointIt++;
>> i++;
>> }
>> cells++;
>>
>> hexElement2->GN = k;
>> m_Solver->el.push_back(itk::fem::FEMP<itk::fem::Element>(hexElement2));
>>
>> break;
>> }
>>
>> case CellType::TRIANGLE_CELL:
>> {
>> if ( Flag ) // To make sure that the nodes are created just once
>> {
>> for (int j=0; j < numofpoints; j++)
>> {
>> m_Input->GetPoint(k, ptr);
>>
>> TriangNodeType::Pointer triangNode;
>> triangNode = new TriangNodeType(pt[0], pt[1], pt[2]);
>> triangNode->GN = j;
>> m_Solver->node.push_back(itk::fem::FEMP<TriangNodeType>(triangNode));
>> }
>> Flag = false;
>> }
>> PointIdIterator pointIt = cellPtr->PointIdsBegin() ;
>> int i=0;
>> triangElement2 = dynamic_cast<TriangElementType*> (triangElement->Clone());
>> while (pointIt != cellPtr->PointIdsEnd() )
>> {
>> triangElement2->SetNode(i, m_Solver->node.Find( *pointIt ));
>> pointIt++; i++;
>> }
>> cells++;
>>
>> triangElement2->GN = k;
>> m_Solver->el.push_back(itk::fem::FEMP<itk::fem::Element>(triangElement2));
>> break;
>> }
>> }
>> }
>>
>> //std::cerr << "Converted Cells " << std::endl;
>>
>> /* #ifdef __sgi
>> // Create the file. This is required on some older sgi's
>> std::ofstream tFile(m_FileName.c_str( ),std::ios::out);
>> tFile.close( );
>> #endif*/
>> /*std::ofstream out;
>> out.open( m_FileName.c_str( ), std::ios::out );
>> m_Solver->Write( out );
>> out.close( );
>> m_Solver->Write(std::cout);*/
>>
>> }
>>}
>>
>>#endif
>>
>>
>>------------------------------------------------------------------------
>>
>>
>>#include <iostream>
>>#include "vtkUnstructuredGridToitkMesh.h"
>>#include <vtkIdList.h>
>>
>>
>>#ifndef vtkDoubleType
>>#define vtkDoubleType double
>>#endif
>>
>>#ifndef vtkFloatingPointType
>># define vtkFloatingPointType vtkFloatingPointType
>>typedef float vtkFloatingPointType;
>>#endif
>>
>>vtkUnstructuredGridToitkMesh
>>::vtkUnstructuredGridToitkMesh()
>>{
>>
>>}
>>
>>
>>vtkUnstructuredGridToitkMesh
>>::~vtkUnstructuredGridToitkMesh()
>>{
>> if (m_Grid)
>> {
>> m_Grid->Delete();
>> }
>>}
>>
>>void
>>vtkUnstructuredGridToitkMesh
>>::SetInput(vtkUnstructuredGrid * grid)
>>{
>> m_Grid = grid;
>> this->ConvertvtkToitk();
>>}
>>
>>vtkUnstructuredGridToitkMesh::MeshType::Pointer
>>vtkUnstructuredGridToitkMesh
>>::GetOutput()
>>{
>> return m_itkMesh;
>>}
>>
>>void
>>vtkUnstructuredGridToitkMesh
>>::ConvertvtkToitk()
>>{
>> m_itkMesh = MeshType::New();
>>
>> // Get the points from vtk
>> vtkPoints* vtkpoints = m_Grid->GetPoints();
>> int numPoints = vtkpoints->GetNumberOfPoints();
>>
>> // Create a compatible point container for the mesh
>> // the mesh is created with a null points container
>> MeshType::PointsContainer::Pointer points =
>> MeshType::PointsContainer::New();
>>
>> // Resize the point container to be able to fit the vtk points
>> points->Reserve(numPoints);
>>
>> // Set the point container on the mesh
>> m_itkMesh->SetPoints(points);
>> for(int i =0; i < numPoints; i++)
>> {
>> double* apoint = vtkpoints->GetPoint(i);
>> float apt[3];
>> apt[0] = static_cast<float>(apoint[0]);
>> apt[1] = static_cast<float>(apoint[1]);
>> apt[2] = static_cast<float>(apoint[2]);
>> m_itkMesh->SetPoint(i, MeshType::PointType(apt));
>> }
>>
>> vtkCellArray* vtkcells = m_Grid->GetCells();
>> MeshType::CellsContainerPointer cells = MeshType::CellsContainer::New();
>>
>>
>> // extract the cell id's from the vtkUnstructuredGrid
>> int numcells = vtkcells->GetNumberOfCells();
>>
>> // extract cell type for each cell
>> int* vtkCellTypes = new int[numcells];
>> int cellId =0;
>> for(; cellId < numcells; cellId++)
>> {
>> vtkCellTypes[cellId] = m_Grid->GetCellType(cellId);
>> }
>>
>> /*
>> std::ofstream out;
>> out.open( "VTKtoITk.txt", std::ios::out );
>> */
>> cells->Reserve(numcells);
>> m_itkMesh->SetCells(cells);
>>
>> vtkIdType npts;
>> vtkIdType* pts;
>> cellId = 0;
>> for(vtkcells->InitTraversal(); vtkcells->GetNextCell(npts, pts); cellId++)
>> {
>> //out << " \n VTK cell ID = " << cellId << " " << vtkCellTypes[cellId] << "\n" << std::endl;
>> MeshType::CellAutoPointer c;
>> switch(vtkCellTypes[cellId])
>> case VTK_HEXAHEDRON:
>> {
>> typedef itk::CellInterface<double, MeshType::CellTraits> CellInterfaceType;
>> typedef itk::HexahedronCell<CellInterfaceType> HexahedronCellType;
>> HexahedronCellType * t = new HexahedronCellType;
>> for ( int i = 0; i < 8; i ++)
>> {
>> //out << "\tVTK Point ids = " << pts[i] << std::endl;
>> t->SetPointId(i, static_cast<unsigned long>(pts[i]));
>> }
>>
>> //t->SetPointIds((unsigned long*)pts);
>> c.TakeOwnership( t );
>> break;
>> }
>> m_itkMesh->SetCell(cellId, c);
>> }
>>
>>}
>>
>>
>>
>>------------------------------------------------------------------------
>>
>>#ifndef __vtkUnstructuredGridToitkMesh_h__
>>
>>#define __vtkUnstructuredGridToitkMesh_h__
>>
>>
>>
>>#include "vtkPoints.h"
>>
>>#include "vtkUnstructuredGrid.h"
>>
>>#include "vtkHexahedron.h"
>>
>>#include "vtkTetra.h"
>>
>>#include "itkDefaultDynamicMeshTraits.h"
>>
>>#include "itkMesh.h"
>>
>>#include "itkAutomaticTopologyMeshSource.h"
>>
>>#include "itkHexahedronCell.h"
>>
>>#include "vtkCellArray.h"
>>
>>
>>
>>/**
>>
>> \class vtkUnstructuredGridToitkMesh
>>
>> \brief
>>
>> \warning
>>
>> \sa
>>
>> */
>>
>>
>>
>>class vtkUnstructuredGridToitkMesh
>>
>>{
>>
>>
>>
>> public:
>>
>>
>>
>> vtkUnstructuredGridToitkMesh( void );
>>
>> ~vtkUnstructuredGridToitkMesh( void );
>>
>>
>>
>> typedef itk::DefaultStaticMeshTraits<double, 3, 3,float,float> MeshTraits;
>>
>> typedef itk::Mesh<double,3, MeshTraits> MeshType;
>>
>>
>>
>> /**
>>
>> Useful Internal Typedefs
>>
>> */
>>
>> typedef itk::AutomaticTopologyMeshSource< MeshType > MeshSourceType;
>>
>> typedef MeshSourceType::Pointer MeshSourceTypePointer;
>>
>>
>>
>>
>>
>> /**
>>
>> The SetInput method provides pointer to the vtkUnstructuredGrid
>>
>> */
>>
>> void SetInput( vtkUnstructuredGrid * grid);
>>
>> MeshType::Pointer GetOutput();
>>
>> void ConvertvtkToitk();
>>
>>
>>
>> MeshType::Pointer m_itkMesh;
>>
>>
>>
>> vtkUnstructuredGrid *m_Grid;
>>
>>
>>
>>
>>
>>};
>>
>>
>>
>>#endif
>>
>>
>>
>>------------------------------------------------------------------------
>>
>>_______________________________________________
>>Insight-users mailing list
>>Insight-users at itk.org
>>http://www.itk.org/mailman/listinfo/insight-users
More information about the Insight-users
mailing list