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

[Francois Budin]

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::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!
> 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/20171009/6528404e/attachment.html>