[ITK-users] itk::Fem::Solver update problem

Dženan Zukić dzenanz at gmail.com
Wed Oct 11 14:06:17 EDT 2017


Hi Yusuf,

as posting attachments on the discourse
<https://discourse.itk.org/t/itk-fem-solver-update-problem/280> forum is
not allowed, here is the updated source code to run on VTK8.

Regards,
Dženan

On Tue, Oct 10, 2017 at 11:36 AM, Yusuf OEZBEK via Insight-users <
insight-users at itk.org> wrote:

> Hello,
>
> Thank you for the answer Francois. Even chaning from raw- to
> smart-pointer for itkMesh did'nt bring any solution. The number of proints,
> cells and cell types are the same for vtkUnstructured grid and itkMesh
> after converting. But I got a new error : "vnl_error_vector_dimension:vnl_matrix::set_row:
> Dimensions [3] and [2] donot match" after the calling m_FemSolver->Update();
> function.
> Do you (or someone else) know, which matrix dimension is meant hier?
>
>
> Same quation is also asked here. But unfortunately without answer:
> https://itk.org/pipermail/insight-users/2009-July/031536.html
>
>
>
> Best regards
>
>
>
> Francois Budin <francois.budin at kitware.com> schrieb am 19:35 Montag,
> 9.Oktober 2017:
>
>
> Hello Yusuf,
>
> I am not sure what the problem is, but here are some ideas you can try to
> understand it or try to solve it:
> 1) Use ITK smart pointers instead of raw pointers. The crash seem to
> indicate that the program was trying to unregister a smart pointer, but
> your mesh is declared as a raw pointer. The issue could come from this.
> 2) Have you tried to look at your data, once converted, to make sure that
> it was converted correctly? Maybe there is a problem during your conversion
> which leads to an empty mesh or something worst.
>
> Hope this helps,
> Francois
>
> PS: As a side note, you may want to join the ITK discourse forum to ask
> your question or help others. The mailing list is still active but more
> people seem to be using the discourse forum now that it is available (
> http://discourse.itk.org/)
>
> On Mon, Oct 9, 2017 at 8:08 AM, Yusuf OEZBEK via Insight-users <
> insight-users at itk.org> wrote:
>
> Dear ITK Users,
>
> According the mail https://public.kitware.com/
> pipermail/insight-users/2007- November/024272.html
> <https://public.kitware.com/pipermail/insight-users/2007-November/024272.html>
> , I have converted my vtkUnstructuredGrid, wich I extracted from a
> vtkPolyData and contains TRIANGLE_CELL's to the itkMesh and than from
> itkMesh to itk::Fem::Object to deform my mesh:
> My problem is: If I call "m_FemSolver->Update();" ,I get a segmentation
> fault, which I really not understand.
>
> Program received signal SIGSEGV, Segmentation fault.
> 0x00007fffb1086dd5 in itk::fem::FEMObject<2u>:: DeepCopy (this=0x4cc5ee0,
>     Copy=0x4cc1780)
>     at /home/toolkits/itk-4.7.2/ Modules/Numerics/FEM/include/
> itkFEMObject.hxx:157
> 157        a->UnRegister();
>
>
> void itkMeshToFemMesh:: updateItkMeshSlot(itk::Mesh< double, 3,
> MeshTraitsType> *itkMesh)
> {
>   typedef itk::fem::FEMObject<3> FemObjectType;
>   FemObjectType::Pointer m_FemObject= FemObjectType::New();
>
>   typedef itk::fem::Solver<3> SolverType;
>   typedef SolverType::Pointer SolverPointerType;
>   SolverPointerType m_FemSolver= SolverType::New();
>
>   // Define the Element Material properties
>   typedef itk::fem:: MaterialLinearElasticity MaterialType;
>   MaterialType::Pointer material= MaterialType::New();
>   material->SetGlobalNumber(0);
>   material->SetYoungsModulus( 200E6);
>   material-> SetCrossSectionalArea(1.0);
>   material->SetThickness(1.0);
>   material->SetMomentOfInertia( 1.0);
>   material->SetPoissonsRatio(0. 2);
>   material-> SetDensityHeatProduct(1.0);
>   m_FemObject->AddNextMaterial( material);
>
>   // Define the element type
>   typedef itk::fem:: Element2DC0LinearTriangularStr ain
> TriangularElementType;
>
>   TriangularElementType::Pointer triangularElement=
> TriangularElementType::New();
>   triangularElement-> SetMaterial(material. GetPointer());
>
>   // Convert mesh points into nodes
>   VectorType point(3);
>   PointType*ptr;
>   PointType pt;
>   ptr= &pt;
>
>   int numOfPoints= itkMesh->GetNumberOfPoints();
>   cout<<"itkMesh numOfPoints: "<<numOfPoints<<endl;
>
>   int numberOfCells= itkMesh->GetNumberOfCells();
>   cout<<"itkMesh numberOfCells: "<<numberOfCells<<endl;
>
>   CellsContainerPointer cellIterator= itkMesh->GetCells();
>   CellIterator cells= cellIterator->Begin();
>
>   bool flag= true;
>
>   for(int k=0; k< numberOfCells; k++)
>   {
>     CellType *cellPtr= cells.Value();
>     cout<<"Cell Value: "<< cells.Value() << " Cell Type= " <<
> cellPtr->GetType()<<endl;
>
>     switch(cellPtr->GetType())
>     {
>       case CellType::TRIANGLE_CELL:
>       {
>         if(flag== true) // To make sure that the nodes are created just
> once
>         {
>           for(int j= 0; j< numOfPoints; j++)
>           {
>             itkMesh->GetPoint(j, ptr);
>
>             typedef TriangularElementType::Node TriangularNodeType;
>             TriangularNodeType::Pointer triangularNode=
> TriangularNodeType::New();
>             point[0]= -1.0;
>             point[1]= 2.0;
>             point[2]= 3.0;
>             triangularNode-> SetCoordinates(point);
>             triangularNode-> SetGlobalNumber(j);
>             m_FemObject->AddNextNode( triangularNode.GetPointer());
>           }
>           flag= false;
>         }
>         PointIdIterator pointIt= cellPtr->PointIdsBegin();
>         int i= 0;
>
>         while(pointIt!= cellPtr->PointIdsEnd())
>         {
>            triangularElement->SetNode(i, m_FemObject->GetNode(*pointIt) );
>           pointIt++;
>           i++;
>         }
>         cells++;
>         triangularElement-> SetGlobalNumber(k);
>         m_FemObject->AddNextElement( triangularElement.GetPointer() );
>         break;
>       }
>     }
>   }
>
>   // Define some Load
>   LoadNodePointerType loadNode =LoadNodeType::New();
>   loadNode->SetElement( triangularElement);
>   loadNode->SetGlobalNumber(3);
>   loadNode->SetNode(1);
>
>   vnl_vector<double> force(2);
>   force[0]= 20.0;
>   force[1]= -20.0;
>   force[2]= 10.0;
>   loadNode->SetForce(force);
>   m_FemObject->AddNextLoad( loadNode);
>
>   //Solve for displacements.
>   m_FemObject->FinalizeMesh();
>   m_FemSolver->SetInput(m_ FemObject);
>   m_FemSolver->Update();     // SEGMENTATION FAULT
>
>   cout<< "Fem Solver Output: "<< m_FemSolver->GetOutput()<< endl;
>
>   const unsigned int invalidId= itk::fem::Element::
> InvalidDegreeOfFreedomID;
>   int numberOfNodes= m_FemSolver->GetInput()-> GetNumberOfNodes();
>
>   for(int i= 0; i< numberOfNodes; i++)
>   {
>     itk::fem::Element::Node:: Pointer node= m_FemSolver->GetInput()->
> GetNode(i);
>     cout<<"Fem Node : "<< node->GetGlobalNumber()<<endl;
>
>     for(unsigned int j= 0, dof; (dof= node->GetDegreeOfFreedom(j))!=
> invalidId; j++)
>     {
>       cout <<"Solver GetSolution: "<< m_FemSolver->GetSolution(dof)< <endl;
>     }
>   }
>
>   // Write the deformed mesh
>   typedef itk::FEMObjectSpatialObject<3> FemObjectSpatialObjectType;
>   FemObjectSpatialObjectType:: Pointer m_FemSpatialObject=
> FemObjectSpatialObjectType:: New();
>   m_FemSpatialObject-> SetFEMObject(m_FemSolver-> GetOutput());
>
>   typedef itk::FEMSpatialObjectWriter<3>
> FemSpatialObjectWriterType;
>   FemSpatialObjectWriterType:: Pointer femSpatialWriter=
> FemSpatialObjectWriterType:: New();
>   femSpatialWriter->SetInput(m_ FemSpatialObject);
>   femSpatialWriter->SetFileName( "deformedMesh.meta");
>   femSpatialWriter->Update();
> }
>
> Thank you for any help!
> The ITK community is transitioning from this mailing list to
> discourse.itk.org. Please join us there!
> ______________________________ __
> Powered by www.kitware.com
>
> Visit other Kitware open-source projects at
> http://www.kitware.com/ opensource/opensource.html
> <http://www.kitware.com/opensource/opensource.html>
>
> Kitware offers ITK Training Courses, for more information visit:
> http://www.kitware.com/ products/protraining.php
> <http://www.kitware.com/products/protraining.php>
>
> Please keep messages on-topic and check the ITK FAQ at:
> http://www.itk.org/Wiki/ITK_ FAQ <http://www.itk.org/Wiki/ITK_FAQ>
>
> Follow this link to subscribe/unsubscribe:
> http://public.kitware.com/ mailman/listinfo/insight-users
> <http://public.kitware.com/mailman/listinfo/insight-users>
>
>
>
>
>
> The ITK community is transitioning from this mailing list to
> discourse.itk.org. Please join us there!
> ________________________________
> Powered by www.kitware.com
>
> Visit other Kitware open-source projects at
> http://www.kitware.com/opensource/opensource.html
>
> Kitware offers ITK Training Courses, for more information visit:
> http://www.kitware.com/products/protraining.php
>
> Please keep messages on-topic and check the ITK FAQ at:
> http://www.itk.org/Wiki/ITK_FAQ
>
> Follow this link to subscribe/unsubscribe:
> http://public.kitware.com/mailman/listinfo/insight-users
>
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://public.kitware.com/pipermail/insight-users/attachments/20171011/637007f4/attachment.html>
-------------- next part --------------
cmake_minimum_required(VERSION 3.8)
 
PROJECT(vtkUGridToFemMesh)

find_package(VTK REQUIRED)
include(${VTK_USE_FILE})

find_package(ITK REQUIRED)
include(${ITK_USE_FILE})


INCLUDE_DIRECTORIES(
  ${vtkCone_SOURCE_DIR}
  ${vtkCone_BINARY_DIR}
)
 

SET(vtkUGridToFemMesh_SRCS
main.cxx
)

add_executable(vtkUGridToFemMesh ${vtkUGridToFemMesh_SRCS})
 
target_link_libraries(vtkUGridToFemMesh ${VTK_LIBRARIES} ${ITK_LIBRARIES})

-------------- next part --------------

#include <iostream>
#include <vtkSmartPointer.h>
#include <vtkPolyDataReader.h>
#include <vtkUnstructuredGrid.h>
#include <vtkDataSetMapper.h>
#include <vtkActor.h>
#include <vtkRenderer.h>
#include <vtkRenderWindow.h>
#include <vtkRenderWindowInteractor.h>
#include <vtkAppendFilter.h>
#include <vtkPoints.h>
#include <vtkHexahedron.h>
#include <vtkTriangle.h>
#include <vtkQuad.h>
#include <vtkTetra.h>
#include <vtkVertex.h>
#include <vtkCellArray.h>

#include <itkDefaultDynamicMeshTraits.h>
#include <itkMesh.h>
#include <itkAutomaticTopologyMeshSource.h>
#include <itkHexahedronCell.h>
#include <itkTriangleCell.h>
#include <itkQuadrilateralCell.h>
#include <itkPolygonCell.h>
#include <itkFEMObject.h>
#include <itkFEMSolver.h>
#include <itkFEMElement3DC0LinearHexahedronStrain.h>
#include <itkFEMElement2DC0LinearTriangularStrain.h>
#include <itkFEMLoadBC.h>
#include <itkFEMObjectSpatialObject.h>
#include <itkFEMSpatialObjectWriter.h>
#include <itkMetaMeshConverter.h>
#include <itkMetaConverterBase.h>
#include <itkMeshSpatialObject.h>


typedef itk::DefaultStaticMeshTraits<double, 3, 3, float, float> MeshTraitsType;
typedef itk::Mesh<double, 3, MeshTraitsType> MeshType;


int main(int, char *[])
{

  /////////////////////////// CONVERT VTK POLY DATA TO UNSTRUCTURED GRID ///////////////////////////
  cout<<"CONVERT VTK POLY DATA TO UNSTRUCTURED GRID"<<endl;

  //Read the content of .ext file
  vtkSmartPointer<vtkPolyDataReader> polyReader= vtkSmartPointer<vtkPolyDataReader>::New();
  polyReader->SetFileName("extracted.ext");
  polyReader->Update();

  vtkSmartPointer<vtkAppendFilter> append= vtkSmartPointer<vtkAppendFilter>::New();
  append->AddInputConnection(polyReader->GetOutputPort());
  append->Update();

  vtkSmartPointer<vtkUnstructuredGrid> uGrid= vtkSmartPointer<vtkUnstructuredGrid>::New();
  uGrid->ShallowCopy(append->GetOutput());

  //Create a mapper and actor
  vtkSmartPointer<vtkDataSetMapper> mapper= vtkSmartPointer<vtkDataSetMapper>::New();
  mapper->SetInputData(uGrid);
  mapper->ScalarVisibilityOff();

  vtkSmartPointer<vtkActor> actor= vtkSmartPointer<vtkActor>::New();
  actor->SetMapper(mapper);

  vtkSmartPointer<vtkRenderer> renderer= vtkSmartPointer<vtkRenderer>::New();
  vtkSmartPointer<vtkRenderWindow> renderWindow= vtkSmartPointer<vtkRenderWindow>::New();
  renderWindow->AddRenderer(renderer);

  vtkSmartPointer<vtkRenderWindowInteractor> renderWindowInteractor= vtkSmartPointer<vtkRenderWindowInteractor>::New();
  renderWindowInteractor->SetRenderWindow(renderWindow);

  renderer->AddActor(actor);
  renderer->SetBackground(0, 0, 0);

  renderWindow->Render();
  renderWindowInteractor->Start();



  /////////////////////////// CONVERT VTK UNSTRUCTURED GRID TO ITK MESH ///////////////////////////
  cout<<"CONVERT VTK UNSTRUCTURED GRID TO ITK MESH"<<endl;
  MeshType::Pointer itkMesh= MeshType::New();

    // Get the points from vtk
    vtkPoints* vtkPoints= uGrid->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
    itkMesh->SetPoints(points);
    for(int i= 0; i< numPoints; i++)
    {
      double *aPoint= vtkPoints->GetPoint(i);
      float itkPoint[3];
      itkPoint[0]= static_cast<float>(aPoint[0]);
      itkPoint[1]= static_cast<float>(aPoint[1]);
      itkPoint[2]= static_cast<float>(aPoint[2]);

      itkMesh->SetPoint(i, MeshType::PointType(itkPoint));
    }

    vtkSmartPointer<vtkCellArray> vtkCells= vtkSmartPointer<vtkCellArray>::New();
    vtkCells= uGrid->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]= uGrid->GetCellType(cellId);
    }

    // Resize the point container to be able to fit the vtk points
    cells->Reserve(numCells);
    itkMesh->SetCells(cells);

    vtkIdType npts;
    vtkIdType *pts;
    cellId= 0;

    for(vtkCells->InitTraversal(); vtkCells->GetNextCell(npts, pts); cellId++)
    {
      //cout<<"VTK cell ID= "<<cellId <<" "<<vtkCellTypes[cellId]<<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++)
  	{
  	  t->SetPointId(i, static_cast<unsigned long>(pts[i]));
  	}
  	  //t->SetPointIds((unsigned long*)pts);
  	  c.TakeOwnership(t);
  	  break;
         }

  	case VTK_TRIANGLE:
  	{
  	  typedef itk::CellInterface<double, MeshType::CellTraits> CellInterfaceType;
  	  typedef itk::TriangleCell<CellInterfaceType> TriangleCellType;

  	  TriangleCellType *t= new TriangleCellType;
  	  for(int i= 0; i< 3; i++)
  	  {
  	    //cout <<"VTK Point ids= "<< pts[i]<<endl;
  	    t->SetPointId(i, static_cast<unsigned long>(pts[i]));
  	  }
  	  c.TakeOwnership(t);
  	  break;
  	}
  	case VTK_POLY_VERTEX:
  	case VTK_LINE:
  	case VTK_POLY_LINE:
  	case VTK_TRIANGLE_STRIP:
  	case VTK_POLYGON:
  	case VTK_PIXEL:
  	case VTK_TETRA:
  	case VTK_VOXEL:
  	case VTK_WEDGE:
  	case VTK_PYRAMID:
  	case VTK_PARAMETRIC_CURVE:
  	case VTK_PARAMETRIC_SURFACE:
  	case VTK_EMPTY_CELL:
  	default:
  	cout<< "Warning unhandled cell type" << vtkCellTypes[cellId]<<endl;
      }
      itkMesh->SetCell(cellId, c);
    }



    /////////////////////////// CONVERT ITK MESH TO FEM MESH ///////////////////////////
    cout<<"CONVERT ITK MESH TO FEM MESH"<<endl;
    typedef itk::fem::FEMObject<3> FemObjectType;
    FemObjectType::Pointer m_FemObject= FemObjectType::New();

     typedef itk::fem::Solver<3> SolverType;
     SolverType::Pointer m_FemSolver= SolverType::New();

     // Define the Element Material properties
     typedef itk::fem::MaterialLinearElasticity MaterialType;
     MaterialType::Pointer material= MaterialType::New();
     material->SetGlobalNumber(0);
     material->SetYoungsModulus(200E6);
     material->SetCrossSectionalArea(1.0);
     material->SetThickness(1.0);
     material->SetMomentOfInertia(1.0);
     material->SetPoissonsRatio(0.2);
     material->SetDensityHeatProduct(1.0);
     m_FemObject->AddNextMaterial(material);

     // Define the element type
     typedef itk::fem::Element2DC0LinearTriangularStrain TriangularElementType;
     TriangularElementType::Pointer triangularElement= TriangularElementType::New();
     triangularElement->SetMaterial(material.GetPointer());

     // Convert mesh points into nodes
     typedef itk::fem::Element::VectorType VectorType;
     VectorType point(3);

     typedef itk::DefaultStaticMeshTraits <double, 3, 3, float, float> MeshTraitsType;
     typedef itk::Mesh <double, 3, MeshTraitsType> MeshType;
     typedef MeshType::PointType PointType;
     PointType*ptr;
     PointType pt;
     ptr= &pt;


     int numOfPoints= itkMesh->GetNumberOfPoints();
     cout<<"itkMesh numOfPoints: "<<numOfPoints<<endl;

     int numberOfCells= itkMesh->GetNumberOfCells();
     cout<<"itkMesh numberOfCells: "<<numberOfCells<<endl;
     /*
     // Access points
     typedef MeshType::PointsContainer::Iterator PointsIterator;
     PointsIterator pointIterator= itkMesh->GetPoints()->Begin();
     PointsIterator endOfPoints= itkMesh->GetPoints()->End();
     MeshType::PointType itkMeshPoints;

     while(pointIterator!= endOfPoints)
     {
       itkMeshPoints= pointIterator.Value();
       cout <<"itkMeshPoints: "<< itkMeshPoints <<endl;
       ++pointIterator;
     }
    */

     for (int j = 0; j< numOfPoints; j++)
     {
         point = itkMesh->GetPoint(j, ptr);

         typedef TriangularElementType::Node TriangularNodeType;
         TriangularNodeType::Pointer triangularNode = TriangularNodeType::New();
         point[0] = -4.0;
         point[1] = 3.0;
         point[2] = 4.0;
         triangularNode->SetCoordinates(point);
         triangularNode->SetGlobalNumber(j);
         m_FemObject->AddNextNode(triangularNode.GetPointer());
         //m_FemObject->RenumberNodeContainer();

         //cout << "j:" << j << " Mesh points: " <<  point[0] << "|" << point[1] << "|" << point[2] << endl;
     }

     typedef MeshType::CellsContainerPointer CellsContainerPointer;
     typedef MeshType::CellsContainer::ConstIterator  CellIterator;
     CellsContainerPointer cellIterator= itkMesh->GetCells();
     CellIterator cellsIter= cellIterator->Begin();

     typedef MeshType::CellType CellType;

     for(int k=0; k< numberOfCells; k++)
     {
       CellType *cellPtr= cellsIter.Value();
       //cout<<"Cell Value: "<< cells.Value() << " Cell Type= " << cellPtr->GetType()<<endl;

       switch(cellPtr->GetType())
       {
        case CellType::TRIANGLE_CELL:
        {
          typedef itk::TriangleCell <CellType> TriangleType;
          typedef TriangleType::PointIdIterator PointIdIterator;
          PointIdIterator pointIt= cellPtr->PointIdsBegin();
          int i= 0;

          while(pointIt!= cellPtr->PointIdsEnd())
          {
            triangularElement->SetNode(i, m_FemObject->GetNode(*pointIt));
            pointIt++;
            i++;
          }
          cellsIter++;

          triangularElement->SetGlobalNumber(k);
          m_FemObject->AddNextElement(triangularElement.GetPointer());
          break;
         }
       }
     }

     // Define boundary condiditon load types
     typedef itk::fem::LoadBC LoadBCType;
     LoadBCType::Pointer loadBc1= LoadBCType::New();
     loadBc1->SetElement(m_FemObject->GetElement(0));
     loadBc1->SetGlobalNumber(0);
     loadBc1->SetDegreeOfFreedom(0);
     loadBc1->SetValue(vnl_vector<double>(1, 0.0));
     m_FemObject->AddNextLoad(loadBc1);

     LoadBCType::Pointer loadBc2= LoadBCType::New();
     loadBc2->SetElement(m_FemObject->GetElement(0));
     loadBc2->SetGlobalNumber(1);
     loadBc2->SetDegreeOfFreedom(1);
     loadBc2->SetValue(vnl_vector<double>(1, 0.0));
     m_FemObject->AddNextLoad(loadBc2);

     LoadBCType::Pointer loadBc3= LoadBCType::New();
     loadBc3->SetElement(m_FemObject->GetElement(1));
     loadBc3->SetGlobalNumber(2);
     loadBc3->SetDegreeOfFreedom(4);
     loadBc3->SetValue(vnl_vector<double>(1, 0.0));
     m_FemObject->AddNextLoad(loadBc3);

     // Define force type
     typedef itk::fem::LoadNode LoadNodeType;
     LoadNodeType::Pointer loadNode= LoadNodeType::New();
     loadNode->SetElement(m_FemObject->GetElement(2));
     loadNode->SetGlobalNumber(3);
     loadNode->SetNode(1);

     vnl_vector<double> force(3);
     force[0]= 20.0;
     force[1]= -20.0;
     force[2]= 10.0;
     loadNode->SetForce(force);
     m_FemObject->AddNextLoad(loadNode);

     //The whole problem is now stored inside the FEM Object class. Now solve for displacements.
     m_FemObject->FinalizeMesh();

     m_FemSolver->SetInput(m_FemObject);
     m_FemSolver->Update();


     cout<< "Fem Solver Output: "<< m_FemSolver->GetOutput()<<endl;

     const unsigned int invalidId= itk::fem::Element::InvalidDegreeOfFreedomID;
     int numberOfNodes= m_FemSolver->GetInput()->GetNumberOfNodes();

     for(int i= 0; i< numberOfNodes; i++)
     {
       itk::fem::Element::Node::Pointer node= m_FemSolver->GetInput()->GetNode(i);
       cout<<"Fem Node : "<< node->GetGlobalNumber()<<endl;

       for(unsigned int j= 0, dof; (dof= node->GetDegreeOfFreedom(j))!= invalidId; j++)
       {
         cout <<"Solve Solution: "<< m_FemSolver->GetSolution(dof)<<endl;
       }
     }

     // Write the deformed mesh
     typedef itk::FEMObjectSpatialObject<3> FemObjectSpatialObjectType;
     FemObjectSpatialObjectType::Pointer m_FemSpatialObject= FemObjectSpatialObjectType::New();
     m_FemSpatialObject->SetFEMObject(m_FemSolver->GetOutput());

     typedef itk::FEMSpatialObjectWriter<3> FemSpatialObjectWriterType;
     FemSpatialObjectWriterType::Pointer femSpatialWriter= FemSpatialObjectWriterType::New();
     femSpatialWriter->SetInput(m_FemSpatialObject);
     femSpatialWriter->SetFileName("deformedMesh.meta");
     femSpatialWriter->Update();

  return EXIT_SUCCESS;
}



More information about the Insight-users mailing list