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

[Yusuf OEZBEK]

Dear ITK Users,

According the mail 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::Element2DC0LinearTriangularStrain 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!