Dear ITK Users,
According the mail 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
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