itk::Fem::Solver update problem

Dear ITK Users,

According the mail , 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 (FEMMesh) to deform my mesh. The number of points, cells and cell types are the same for vtkUnstructured grid and itkMesh after converting.

My problem is: If I call “m_FemSolver->Update();” ,I get a segmentation fault, which I really not understand:

vnl_error_vector_dimension:vnl_matrix::set_row: Dimensions [3] and [2] do not match.

which matrix dimension is meant hier?

Same quation is also asked here. But unfortunately without answer:

My code looks like:

void itkMeshToFemMesh:: convertItkMeshToFemMesh(itk::Mesh< double, 3, MeshTraitsType>::Pointer 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->SetYoungsModulus( 200E6);
  material-> SetCrossSectionalArea(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 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;

      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) );
        triangularElement-> SetGlobalNumber(k);
        m_FemObject->AddNextElement( triangularElement.GetPointer() );

 // Define boundary condiditon load types
 LoadBCType::Pointer loadBc1= LoadBCType::New();
 loadBc1->SetValue(vnl_vector<double>(1, 0.0));

 LoadBCType::Pointer loadBc2= LoadBCType::New();
 loadBc2->SetValue(vnl_vector<double>(1, 0.0));

 LoadBCType::Pointer loadBc3= LoadBCType::New();
 loadBc3->SetValue(vnl_vector<double>(1, 0.0));

 // Define force type
 LoadNodePointerType loadNode= LoadNodeType::New();

 vnl_vector<double> force(2);
 force[0]= 20.0;
 force[1]= -20.0;

 //Solve for displacements.
 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");

Thanks for any help

This is two-dimensional, and might need to be 3D.

If changing this does not help, can you provide a runnable example with main function and an input file?

Changing the force dimension did not help. Under the link is an example and input file:
To compile it, I use VTK 5.10.1 ,ITK 4.7.2 and Cmake

Thank you.

In file c:\Libs\ITK-4.12.0\Modules\Numerics\FEM\src\itkFEMElementBase.cxx, on line 277 const unsigned int Ndims = this->GetNumberOfSpatialDimensions(); Ndims is 2. That is correct, as all the elements are triangles. A few lines below (line 283) exception occurs, as Nn is 3.

This probably means that you can’t run 3D FEM with 2D elements (triangles). I think that you need 3D elements: tetrahedrons, hexahedrons, cubes.

Hi Dženan

Thank you for the information. I have converted my vtkUnstructuredGrid, which was a VTK_TRIANGLE object to the VTK_TETRA object using vtkDelunay3D class. The output is again a vtkUnstructuredGrid but as TETRAHEDRON. See “extracted2.ext”. I have a also used VTK_TETRAHEDRON cell type in itk to convert it from itkMesh to itkFemMesh and used itk::fem::Element3DC0LinearTetrahedronStrain instead of itk::fem::Element2DC0LinearTriangularStrain. See “main.cxx” Unfortunately I get now error:

/itk-4.7.2/Modules/ThirdParty/VNL/src/vxl/core/vnl/algo/vnl_qr.txx: vnl_qr<T>::solve() : matrix is rank-deficient by 3 

which occurs in:

at /toolkits/itk-4.7.2/Modules/Numerics/FEM/include/itkFEMObject.hxx:157
157	    a->UnRegister();

With this source code, I still get “vnl_error_vector_dimension:vnl_matrix::set_row: Dimensions [3] and [2] do not match.” in the same place as before. Are you sure you provided the right source code?

If I debug with gdb ./vtkUGridToFemMesh
I get:

/home/csad6517/toolkits/itk-4.7.2/Modules/ThirdParty/VNL/src/vxl/core/vnl/algo/vnl_qr.txx: vnl_qr<T>::solve() : matrix is rank-deficient by 3

Program received signal SIGSEGV, Segmentation fault.
0x000000000055f2c2 in itk::fem::FEMObject<3u>::DeepCopy (this=0x1ba33b0,  Copy=0x1b9eaa0)
at /home/csad6517/toolkits/itk-4.7.2/Modules/Numerics/FEM/include/itkFEMObject.hxx:157
157        a->UnRegister();

(gdb) bt
#0  0x000000000055f2c2 in itk::fem::FEMObject<3u>::DeepCopy (this=0x1ba33b0, Copy=0x1b9eaa0) at /home/csad6517/toolkits/itk-4.7.2/Modules/Numerics/FEM/include/itkFEMObject.hxx:157
#1  0x0000000000547105 in itk::fem::Solver<3u>::RunSolver (this=0x1b9ed70) at /home/csad6517/toolkits/itk-4.7.2/Modules/Numerics/FEM/include/itkFEMSolver.hxx:624
#2  0x0000000000545d1f in itk::fem::Solver<3u>::GenerateData (this=0x1b9ed70) at /home/csad6517/toolkits/itk-4.7.2/Modules/Numerics/FEM/include/itkFEMSolver.hxx:178
#3  0x00007ffff14536cd in itk::ProcessObject::UpdateOutputData(itk::DataObject*) () from /home/csad6517    /toolkits/itk-4.7.2_bin/lib/
#4  0x00007ffff146a135 in itk::DataObject::Update() ()from /home/csad6517/toolkits/itk-4.7.2_bin   /lib/
#5  0x0000000000500d9c in main () at /home/csad6517/workspace/vtkUGridToFemMesh/main.cxx:329

Can you download it again please:

Compiler: : gcc (Debian 4.9.2-10) 4.9.2
VTK: 5.10.1
ITK: 4.7.2

In file c:\Libs\ITK-4.12.0\Modules\Numerics\FEM\include\itkFEMObject.hxx, on line 156 a = ObjectFactoryBase::CreateInstance( elCopy->GetNameOfClass() );, a ends up being NULL. Of course, a->UnRegister(); on the next line causes the crash.

I would have expected that class Element gets created by the ObjectFactoryBase. Maybe this part of the code got broken due to some code changes in ObjectFactoryBase, LightObject or something similar. @matt.mccormick @blowekamp @hjmjohnson Can you provide more information?

@jhlegarreta Is function template <unsigned int VDimension> void FEMObject<VDimension>::DeepCopy( FEMObject *Copy) covered by regression testing?

It is recommended to update ITK to the latest release before spending further effort.

This also happens with a recent git version (which I have right now).

I have also tried another cell type (HEXAHEDRON) with itk::fem::Element3DC0LinearHexahedronStrain with newest itk version and got same issue.

Has anyone really no idea about this issue?

I debugged you example, and itkObjectFactoryBase.cxx did not have any factories registered, therefore could not make the element required in itkFEMObject.hxx on line 156. Taking inspiration from itkFEMFactoryTest.cxx, the solution to crash on line 157 is adding the following code in your main.cxx:

#include "itkFEMFactory.h" //add before main
itk::FEMFactoryBase::GetFactory()->RegisterDefaultTypes(); //before Update();

Was ist possible for you to get any solution? I get only 0 from m_FemSolver->GetSolution(…)

I got zeroes too, but that is not surprising given the warnings about mesh quality from VTK and “matrix is rank-deficient” by vnl_qr<T>::solve.

have I this warnings “/itk-4.7.2/Modules/ThirdParty/VNL/src/vxl/core/vnl/algo/vnl_qr.txx: vnl_qr::solve() : matrix is rank-deficient by 3” because my mesh contains few cells therefore few nodes in tetrahedrons?

I don’t have much experience with FEM simulation. Now that the crash is solved, you should start a new topic about higher level issues with your example, to give a better chance to other people to answer.