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

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})

#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;
}