itk FEM Solver

Algorithms
,
1

Dear ITK community,

I want to deform a 3D model using itk’s FEM module. For this I have implemented the following workflow:

  1. Read CT data-set and create 3D iso-surface using vtkMarchingCubes (output is polydata)
  2. Extracted the interested surface using vtkHardwareSelector (output is 2D polydata)
  3. Converted the extracted surface to the unstructured gird using vtkUnstructuredGrid and vtkDelaunay3D in order to get a 3D tetrahedron as unstructured grid (output is 3D unstructured grid)
  4. Converted the unstructured grid to the itk mesh, which contains tetrahedron cells (output is itkMesh object)
  5. To convert the itk mesh to the FEM model, I have created FEM object (FEMObject<3>), which contains a material (MaterialLinearElasticity), elements and nodes (Element3DC0LinearTetrahedronStrain), boundary conditions within force (LoadBC) and to generate a solution for the formulation, I have created a solver (FEMSolver<3>).

My problem is:
If I call the update of solver “solver->Update()” I get the 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” and generated solutions are just zeros.

Thank you for any help

1 Like
2

Ok I found my problem. It was by creating the nodes. For my case; the node coordinates should be defined from itkMesh and saved in a 3D Vector. Then set as coordinate for nodes.

 typedef MeshType::PointType PointType;
 typedef itk::fem::Element::VectorType VectorType;
 VectorType vector3d(3);
 PointType point3d(3);
 for(int j= 0; j< numOfPoints; j++)
 {
   vector3d= itkMesh->GetPoint(j, &point3d);

  TetrahedronNodeType::Pointer tetrahedronNode= TetrahedronNodeType::New();
  vector3d[0]= point3d[0];
  vector3d[1]= point3d[1];
  vector3d[2]= point3d[2];
  
  tetrahedronNode->SetCoordinates(vector3d);
  tetrahedronNode->SetGlobalNumber(j);
  m_FemObject->AddNextNode(tetrahedronNode.GetPointer());
}

Now I want to ask. Is there any possibility to convert the deformed fem object to the vtk (e.g. polydata or unstructured grid), without saving it as meta object?

3

Thanks for sharing your solution.

Of course you could reverse the conversion process. This example should be helpful.

4

I have another question. If I want to replace the old node coordinates with the new one after solution, in order to visualize the deformation,I get zeros as solution values in (m_ls->GetSolutionValue(node->GetDegreeOfFreedom(j)).

The code-block:

//m_Fem3dSolver->UpdateDisplacements();  
int numberOfNodes= m_Fem3dSolver->GetInput()->GetNumberOfNodes();
const unsigned int invalidId= itk::fem::Element::InvalidDegreeOfFreedomID;
for(int i= 0; i< numberOfNodes; i++)
{
  TetrahedronNodeType::Pointer node= m_Fem3dSolver->GetInput()->GetNode(i);
  std::cout<<"TETRAHEDRON_CELL Fem Node : "<< node->GetGlobalNumber() << " Coordinates: "<< node-
   >GetCoordinates()<<std::endl;

   for(unsigned int j= 0, dof; (dof= node->GetDegreeOfFreedom(j))!= invalidId; j++)
   {
     std::cout <<"TETRAHEDRON_CELL Solver Solution: "<< m_Fem3dSolver->GetSolution(dof)<<std::endl;
   }
 }

gives me following output, which shows I have some replacements for node 1 and 2:

TETRAHEDRON_CELL Fem Node : 0 Coordinates: -60.090999603271484375 -139.6909942626953125 209.5
TETRAHEDRON_CELL Fem Node : 1 Coordinates: -59.80970001220703125 -140.121002197265625 209.5
TETRAHEDRON_CELL Solver Solution: 0
TETRAHEDRON_CELL Solver Solution: 0
TETRAHEDRON_CELL Solver Solution: 0.000790482550413578008603743274563
TETRAHEDRON_CELL Fem Node : 2 Coordinates: -59.785198211669921875 -139.6909942626953125 '
208.7949981689453125
TETRAHEDRON_CELL Solver Solution: -0.000215714426963713613166762073092
TETRAHEDRON_CELL Solver Solution: -9.37249422155208276354054763857e-05
TETRAHEDRON_CELL Solver Solution: 0.000718805934793648290855039295622
TETRAHEDRON_CELL Fem Node : 3 Coordinates: -60.214801788330078125 -139.261993408203125 
209.11700439453125
TETRAHEDRON_CELL Fem Node : 4 Coordinates: -60.36569976806640625 -139.261993408203125 209.5

if I want to update the displacements like in (itkFEMSolver.hxx 610-651):

//itk::fem::LinearSystemWrapper::Pointer wrapper;
VectorType replacedNodeVector(3);

for(int i= 0; i< numberOfNodes; i++)
{
  TetrahedronNodeType::Pointer node= m_Fem3dSolver->GetInput()->GetNode(i);
 //itk::fem::Element::Node::Pointer node= m_Fem3dSolver->GetInput()->GetNode(i);

   for(unsigned int j= 0; j< 3; j++)
  {
    std::cout<<"TETRAHEDRON_CELL Fem Solution Vector: " << m_ls->GetSolutionValue(node->GetDegreeOfFreedom(j))<<std::endl;
    replacedNodeVector[j] = node->GetCoordinates()[j] + m_ls->GetSolutionValue(node->GetDegreeOfFreedom(j));
  }
  node->SetCoordinates(replacedNodeVector);
  std::cout<<"TETRAHEDRON_CELL Fem New Node : "<< node->GetGlobalNumber() << " Coordinates: "<< node->GetCoordinates()<<std::endl;
}

I get as output:

TETRAHEDRON_CELL Fem Solution Vector: 0
TETRAHEDRON_CELL Fem Solution Vector: 0
TETRAHEDRON_CELL Fem Solution Vector: 0
TETRAHEDRON_CELL Fem New Node : 0 Coordinates: -60.090999603271484375 -139.6909942626953125 209.5
TETRAHEDRON_CELL Fem Solution Vector: 0
TETRAHEDRON_CELL Fem Solution Vector: 0
TETRAHEDRON_CELL Fem Solution Vector: 0
TETRAHEDRON_CELL Fem New Node : 1 Coordinates: -59.80970001220703125 -140.121002197265625 209.5
TETRAHEDRON_CELL Fem Solution Vector: 0
TETRAHEDRON_CELL Fem Solution Vector: 0
TETRAHEDRON_CELL Fem Solution Vector: 0
TETRAHEDRON_CELL Fem New Node : 2 Coordinates: -59.785198211669921875 -139.6909942626953125 208.7949981689453125
TETRAHEDRON_CELL Fem Solution Vector: 0
TETRAHEDRON_CELL Fem Solution Vector: 0
TETRAHEDRON_CELL Fem Solution Vector: 0
TETRAHEDRON_CELL Fem New Node : 3 Coordinates: -60.214801788330078125 -139.261993408203125 209.11700439453125
TETRAHEDRON_CELL Fem Solution Vector: 0
TETRAHEDRON_CELL Fem Solution Vector: 0
TETRAHEDRON_CELL Fem Solution Vector: 0
TETRAHEDRON_CELL Fem New Node : 4 Coordinates: -60.36569976806640625 -139.261993408203125 209.5

As you see the old and new node coordinates are the same, because the solution vector gives me zeros, although I have some displacements in solution.

Or should I may be fill the replacedNodeVector with the old node coordinate + Solution for the specified node. like:

replacedNodeVector[j] = node->GetCoordinates()[j] +  m_Fem3dSolver->GetSolution(node->GetDegreeOfFreedom(j));

So I get now different coordinates for node 1 and 2

Where can be my mistake?

Thank you.

5

If you want to visualize mesh deformation, why do you want to change it?

To have a chance of somebody debug this for you, you need to provide a runnable example with input data. Also, make it a new thread for better visibility on the forum.