[ITK-dev] itk::fem updating displacements

[Yusuf OEZBEK]

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:vtkUGridToFemMesh.zip
  
|  
|   
|   
|   |    |

   |

  |
|  
|    |  
vtkUGridToFemMesh.zip
 Shared with Dropbox  |   |

  |

  |

 
Thank you for any help!
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://public.kitware.com/pipermail/insight-developers/attachments/20171109/9820a208/attachment.html>