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

Tue Oct 10 11:36:48 EDT 2017

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 , 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

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/20171010/378f272e/attachment.html>