[ITK-users] itk::Fem::Solver update problem

Dženan Zukić

2017-10-11 18:06:17 UTC

Hi Yusuf,

as posting attachments on the discourse
<https://discourse.itk.org/t/itk-fem-solver-update-problem/280> forum is
not allowed, here is the updated source code to run on VTK8.

Regards,
DÅŸenan

On Tue, Oct 10, 2017 at 11:36 AM, Yusuf OEZBEK via Insight-users <

Hello,
Thank you for the answer Francois. Even chaning from raw- to
smart-pointer for itkMesh did'nt bring any solution. The number of proints,
cells and cell types are the same for vtkUnstructured grid and itkMesh
Dimensions [3] and [2] donot match" after the calling m_FemSolver->Update();
function.
Do you (or someone else) know, which matrix dimension is meant hier?
https://itk.org/pipermail/insight-users/2009-July/031536.html
Best regards
Hello Yusuf,
I am not sure what the problem is, but here are some ideas you can try to
1) Use ITK smart pointers instead of raw pointers. The crash seem to
indicate that the program was trying to unregister a smart pointer, but
your mesh is declared as a raw pointer. The issue could come from this.
2) Have you tried to look at your data, once converted, to make sure that
it was converted correctly? Maybe there is a problem during your conversion
which leads to an empty mesh or something worst.
Hope this helps,
Francois
PS: As a side note, you may want to join the ITK discourse forum to ask
your question or help others. The mailing list is still active but more
people seem to be using the discourse forum now that it is available (
http://discourse.itk.org/)
On Mon, Oct 9, 2017 at 8:08 AM, Yusuf OEZBEK via Insight-users <
Dear ITK Users,
According the mail https://public.kitware.com/
pipermail/insight-users/2007- November/024272.html
<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
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:: 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())
{
{
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;
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
<http://www.kitware.com/opensource/opensource.html>
http://www.kitware.com/ products/protraining.php
<http://www.kitware.com/products/protraining.php>
http://www.itk.org/Wiki/ITK_ FAQ <http://www.itk.org/Wiki/ITK_FAQ>
http://public.kitware.com/ mailman/listinfo/insight-users
<http://public.kitware.com/mailman/listinfo/insight-users>
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
http://www.kitware.com/products/protraining.php
http://www.itk.org/Wiki/ITK_FAQ
http://public.kitware.com/mailman/listinfo/insight-users