itk::fem updating displacements

Algorithms
1

Dear ITK users,

I have another question. In order to visualize the deformation, I want to replace the old node coordinates with the new one after generating the solution, like in UpdateDisplacements() function in itkFemSolver.hxx 610-651
The code-block, which generates solution:

m_FemSolver->Update();
//m_FemSolver->UpdateDisplacements();
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);
  std::cout<<"FEM Mesh Node: "<< node->GetGlobalNumber() << " Coordinates: "<< node->GetCoordinates()<<std::endl;
  
  for(unsigned int j= 0, dof; (dof= node->GetDegreeOfFreedom(j))!= invalidId; j++)
  {
    cout <<"FEM Mesh Solution: "<< m_FemSolver->GetSolution(dof)<<endl;
  }
}

gives me following output, that shows, I have some solutions for nodes 69 and 70:

FEM Mesh Node: 68 Coordinates: -61.1098 -139.262 211.9
FEM Mesh Node: 69 Coordinates: -61.0742 -139.319 211.9
FEM Mesh Solution: 82.6806
FEM Mesh Solution: -32.5517
FEM Mesh Solution: 6.36021
FEM Mesh Node: 70 Coordinates: -61.3429 -138.832 211.9
FEM Mesh Solution: 0.798099
FEM Mesh Solution: 37.6971
FEM Mesh Solution: -4.36149
FEM Mesh Node: 71 Coordinates: -61.5039 -138.402 211.509
FEM Mesh Node: 72 Coordinates: -61.6232 -138.402 211.9

if I want to update the displacements:

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

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

I get as output:

FEM Mesh New Node: 67 Coordinates: -61.0742 -139.262 211.797 
FEM Mesh Solution Vector: 0
FEM Mesh Solution Vector: 0
FEM Mesh Solution Vector: 0
FEM Mesh New Node: 68 Coordinates: -61.1098 -139.262 211.9
FEM Mesh Solution Vector: 82.6806
FEM Mesh Solution Vector: -32.5517
FEM Mesh Solution Vector: 6.36021
FEM Mesh New Node: 69 Coordinates: 21.6064 -171.871 218.26
FEM Mesh Solution Vector: 0.798099
FEM Mesh Solution Vector: 37.6971
FEM Mesh Solution Vector: -4.36149
FEM Mesh New Node: 70 Coordinates: -60.5448 -101.135 207.539
FEM Mesh Solution Vector: 0
FEM Mesh Solution Vector: 0
FEM Mesh Solution Vector: 0

My questions are:
Do I do it correct for the replacements, if I sum the old node coordinate with the generated solution vector for that spacial node?:

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

In my test application, I apply some force at node 1( loadNode->SetNode(1); ) but in generated solution it shows that the force doesn’t act on node 1 and nodes around it but for nodes 68, 69, 82, and 83. All the other solutions are zero. My test mesh contains 86 nodes and 162 elements. If I want to apply the force except node 1,eg at node 20, I get then segmentation fault.

typedef itk::fem::LoadNode LoadNodeType;
LoadNodeType::Pointer loadNode= LoadNodeType::New();
loadNode->SetElement(tetrahedronElement.GetPointer());
loadNode->SetGlobalNumber(3);
loadNode->SetNode(1);

vnl_vector<double> force(3);
force[0]= 110.0;
force[1]= 20.0;
force[2]= 31.0;
loadNode->SetForce(force);
m_FemObject->AddNextLoad(loadNode);

For detailed information please see runnable test code:


Thank you.