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 to deform my mesh:My problem is: If I call
"m_FemSolver->Update();" ,I get a segmentation fault, which I really not
understand.
Program received signal SIGSEGV, Segmentation fault.
0x00007fffb1086dd5 in itk::fem::FEMObject<2u>::DeepCopy (this=0x4cc5ee0,
Copy=0x4cc1780)
at
/home/toolkits/itk-4.7.2/Modules/Numerics/FEM/include/itkFEMObject.hxx:157
157 a->UnRegister();
void itkMeshToFemMesh::updateItkMeshSlot(itk::Mesh<double, 3, MeshTraitsType>
*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::Element2DC0LinearTriangularStrain 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 some Load
LoadNodePointerType loadNode =LoadNodeType::New();
loadNode->SetElement(triangularElement);
loadNode->SetGlobalNumber(3);
loadNode->SetNode(1);
vnl_vector<double> force(2);
force[0]= 20.0;
force[1]= -20.0; force[2]= 10.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();
}
Thank you for any help!
The ITK community is transitioning from this mailing list to discourse.itk.org.
Please join us there!
______________________________________
Powered by www.kitware.com
Visit other Kitware open-source projects at
http://www.kitware.com/opensource/opensource.html
Kitware offers ITK Training Courses, for more information visit:
http://kitware.com/products/protraining.php
Please keep messages on-topic and check the ITK FAQ at:
http://www.itk.org/Wiki/ITK_FAQ
Follow this link to subscribe/unsubscribe:
http://public.kitware.com/mailman/listinfo/insight-developers
