# itk::Fem::Solver update problem

**URL:** <https://discourse.itk.org/t/itk-fem-solver-update-problem/280>\
**Category:** Community\
**Created:** [October 11, 2017, 8:46am UTC](https://discourse.itk.org/t/itk-fem-solver-update-problem/280 "2017-10-11T08:46:20Z")\
**Posts on this page:** 17\
**Page:** 1

<div class="post-metadata">

**Author:** ![nasil122002](https://discourse.itk.org/letter_avatar_proxy/v4/letter/n/48db29/32.png) [@nasil122002](https://discourse.itk.org/u/nasil122002)\
**Post date:** [October 11, 2017, 8:46am UTC](https://discourse.itk.org/t/itk-fem-solver-update-problem/280/1 "2017-10-11T08:46:23Z")

</div>

Dear ITK Users,

According the mail [https://public.kitware.com/pipermail/insight-users/2007-November/024272.html](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 (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: [https://itk.org/pipermail/insight-users/2009-July/031536.html](https://itk.org/pipermail/insight-users/2009-July/031536.html)

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->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 boundary condiditon load types
 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
 LoadNodePointerType loadNode= LoadNodeType::New();
 loadNode->SetElement(m_FemObject->GetElement(2));
 loadNode->SetGlobalNumber(3);
 loadNode->SetNode(1);

 vnl_vector<double> force(2);
 force[0]= 20.0;
 force[1]= -20.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();
}

```

Thanks for any help

---

<div class="post-metadata">

**Author:** ![dzenanz](https://discourse.itk.org/user_avatar/discourse.itk.org/dzenanz/32/1093_2.png) [@dzenanz](https://discourse.itk.org/u/dzenanz)\
**Post date:** [October 11, 2017, 1:49pm UTC](https://discourse.itk.org/t/itk-fem-solver-update-problem/280/2 "2017-10-11T13:49:24Z")

</div>

> [@nasil122002](#):
>
> vnl\_vector\<double\> force(2);

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?

---

<div class="post-metadata">

**Author:** ![nasil122002](https://discourse.itk.org/letter_avatar_proxy/v4/letter/n/48db29/32.png) [@nasil122002](https://discourse.itk.org/u/nasil122002)\
**Post date:** [October 11, 2017, 4:01pm UTC](https://discourse.itk.org/t/itk-fem-solver-update-problem/280/3 "2017-10-11T16:01:38Z")

</div>

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

> **[Dropbox - Link not found](https://www.dropbox.com/s/xqvfifmy49q3ywy/vtkUGridToFemMesh.zip?dl=0)**
>
> Dropbox is a free service that lets you bring your photos, docs, and videos anywhere and share them easily. Never email yourself a file again!

Thank you.

---

<div class="post-metadata">

**Author:** ![dzenanz](https://discourse.itk.org/user_avatar/discourse.itk.org/dzenanz/32/1093_2.png) [@dzenanz](https://discourse.itk.org/u/dzenanz)\
**Post date:** [October 11, 2017, 5:49pm UTC](https://discourse.itk.org/t/itk-fem-solver-update-problem/280/4 "2017-10-11T17:49:58Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![nasil122002](https://discourse.itk.org/letter_avatar_proxy/v4/letter/n/48db29/32.png) [@nasil122002](https://discourse.itk.org/u/nasil122002)\
**Post date:** [October 12, 2017, 1:25pm UTC](https://discourse.itk.org/t/itk-fem-solver-update-problem/280/5 "2017-10-12T13:25:31Z")

</div>

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();

```

> **[Dropbox - Link not found](https://www.dropbox.com/s/g72x6vcnhre1pw0/vtkUGridToFemMesh.zip?dl=0)**
>
> Dropbox is a free service that lets you bring your photos, docs, and videos anywhere and share them easily. Never email yourself a file again!

---

<div class="post-metadata">

**Author:** ![dzenanz](https://discourse.itk.org/user_avatar/discourse.itk.org/dzenanz/32/1093_2.png) [@dzenanz](https://discourse.itk.org/u/dzenanz)\
**Post date:** [October 12, 2017, 2:24pm UTC](https://discourse.itk.org/t/itk-fem-solver-update-problem/280/6 "2017-10-12T14:24:47Z")

</div>

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?

---

<div class="post-metadata">

**Author:** ![nasil122002](https://discourse.itk.org/letter_avatar_proxy/v4/letter/n/48db29/32.png) [@nasil122002](https://discourse.itk.org/u/nasil122002)\
**Post date:** [October 12, 2017, 3:29pm UTC](https://discourse.itk.org/t/itk-fem-solver-update-problem/280/7 "2017-10-12T15:29:06Z")

</div>

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/libITKCommon-4.7.so.1
#4 0x00007ffff146a135 in itk::DataObject::Update() ()from /home/csad6517/toolkits/itk-4.7.2_bin /lib/libITKCommon-4.7.so.1
#5 0x0000000000500d9c in main () at /home/csad6517/workspace/vtkUGridToFemMesh/main.cxx:329
(gdb) 

```

Can you download it again please:

> **[Dropbox - Link not found](https://www.dropbox.com/s/xi64zsqlnnvt6b6/vtkUGridToFemMesh.zip?dl=0)**
>
> Dropbox is a free service that lets you bring your photos, docs, and videos anywhere and share them easily. Never email yourself a file again!

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

---

<div class="post-metadata">

**Author:** ![dzenanz](https://discourse.itk.org/user_avatar/discourse.itk.org/dzenanz/32/1093_2.png) [@dzenanz](https://discourse.itk.org/u/dzenanz)\
**Post date:** [October 12, 2017, 10:04pm UTC](https://discourse.itk.org/t/itk-fem-solver-update-problem/280/8 "2017-10-12T22:04:29Z")

</div>

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?

---

<div class="post-metadata">

**Author:** ![matt.mccormick](https://discourse.itk.org/user_avatar/discourse.itk.org/matt.mccormick/32/7_2.png) [@matt.mccormick](https://discourse.itk.org/u/matt.mccormick)\
**Post date:** [October 12, 2017, 10:15pm UTC](https://discourse.itk.org/t/itk-fem-solver-update-problem/280/9 "2017-10-12T22:15:09Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![dzenanz](https://discourse.itk.org/user_avatar/discourse.itk.org/dzenanz/32/1093_2.png) [@dzenanz](https://discourse.itk.org/u/dzenanz)\
**Post date:** [October 12, 2017, 10:29pm UTC](https://discourse.itk.org/t/itk-fem-solver-update-problem/280/10 "2017-10-12T22:29:12Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![nasil122002](https://discourse.itk.org/letter_avatar_proxy/v4/letter/n/48db29/32.png) [@nasil122002](https://discourse.itk.org/u/nasil122002)\
**Post date:** [October 13, 2017, 12:08pm UTC](https://discourse.itk.org/t/itk-fem-solver-update-problem/280/11 "2017-10-13T12:08:31Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![nasil122002](https://discourse.itk.org/letter_avatar_proxy/v4/letter/n/48db29/32.png) [@nasil122002](https://discourse.itk.org/u/nasil122002)\
**Post date:** [October 17, 2017, 2:41pm UTC](https://discourse.itk.org/t/itk-fem-solver-update-problem/280/12 "2017-10-17T14:41:52Z")

</div>

Has anyone really no idea about this issue?

---

<div class="post-metadata">

**Author:** ![dzenanz](https://discourse.itk.org/user_avatar/discourse.itk.org/dzenanz/32/1093_2.png) [@dzenanz](https://discourse.itk.org/u/dzenanz)\
**Post date:** [October 17, 2017, 7:01pm UTC](https://discourse.itk.org/t/itk-fem-solver-update-problem/280/13 "2017-10-17T19:01:56Z")

</div>

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](https://github.com/InsightSoftwareConsortium/ITK/blob/master/Modules/Numerics/FEM/test/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();
```

---

<div class="post-metadata">

**Author:** ![nasil122002](https://discourse.itk.org/letter_avatar_proxy/v4/letter/n/48db29/32.png) [@nasil122002](https://discourse.itk.org/u/nasil122002)\
**Post date:** [October 18, 2017, 3:04pm UTC](https://discourse.itk.org/t/itk-fem-solver-update-problem/280/14 "2017-10-18T15:04:10Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![dzenanz](https://discourse.itk.org/user_avatar/discourse.itk.org/dzenanz/32/1093_2.png) [@dzenanz](https://discourse.itk.org/u/dzenanz)\
**Post date:** [October 18, 2017, 4:59pm UTC](https://discourse.itk.org/t/itk-fem-solver-update-problem/280/15 "2017-10-18T16:59:31Z")

</div>

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`.

---

<div class="post-metadata">

**Author:** ![nasil122002](https://discourse.itk.org/letter_avatar_proxy/v4/letter/n/48db29/32.png) [@nasil122002](https://discourse.itk.org/u/nasil122002)\
**Post date:** [October 18, 2017, 5:15pm UTC](https://discourse.itk.org/t/itk-fem-solver-update-problem/280/16 "2017-10-18T17:15:27Z")

</div>

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?

---

<div class="post-metadata">

**Author:** ![dzenanz](https://discourse.itk.org/user_avatar/discourse.itk.org/dzenanz/32/1093_2.png) [@dzenanz](https://discourse.itk.org/u/dzenanz)\
**Post date:** [October 18, 2017, 5:43pm UTC](https://discourse.itk.org/t/itk-fem-solver-update-problem/280/17 "2017-10-18T17:43:28Z")

</div>

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.