Skip to content

[SofaMiscFem] Remove TopologyHandler in FEM to use TopologyData callbacks (part 4) #2390

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
fredroy merged 7 commits into from
Oct 15, 2021
Merged

[SofaMiscFem] Remove TopologyHandler in FEM to use TopologyData callbacks (part 4) #2390

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
7 commits
Select commit Hold shift + click to select a range
3203da4
[SofaMiscFem] Remove TopologyHandler in FEM to directly use TopologyD…
epernod Jun 8, 2021
9cf658a
[SofaGeneralSimpleFem] Remove TopologyHandler in TetrahedralTensorMas…
epernod Oct 8, 2021
a78c394
Fix gcc compilation
epernod Oct 11, 2021
2e2bfff
Fix gcc compilation
epernod Oct 11, 2021
1a4e92e
Fix compilation
epernod Oct 11, 2021
ab492ce
[SofaMiscFem] Fix TetrahedralTensorMassForceField init. The tetrahedr…
epernod Oct 12, 2021
69f4c36
remove commented code
fredroy Oct 15, 2021
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
34 changes: 7 additions & 27 deletions modules/SofaMiscFem/src/SofaMiscFem/FastTetrahedralCorotationalForceField.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -105,35 +105,17 @@ class FastTetrahedralCorotationalForceField : public core::behavior::ForceField<
}
};

class FTCFTetrahedronHandler : public topology::TopologyDataHandler<core::topology::BaseMeshTopology::Tetrahedron, sofa::type::vector<TetrahedronRestInformation> >
{
public:
typedef typename FastTetrahedralCorotationalForceField<DataTypes>::TetrahedronRestInformation TetrahedronRestInformation;

using Index = sofa::Index;

FTCFTetrahedronHandler(FastTetrahedralCorotationalForceField<DataTypes>* ff,
topology::TetrahedronData<sofa::type::vector<TetrahedronRestInformation> >* data )
:topology::TopologyDataHandler<core::topology::BaseMeshTopology::Tetrahedron, sofa::type::vector<TetrahedronRestInformation> >(data)
,ff(ff)
{

}

void applyCreateFunction(Index, TetrahedronRestInformation &t,
const core::topology::BaseMeshTopology::Tetrahedron&,
const sofa::type::vector<Index> &,
const sofa::type::vector<double> &);

protected:
FastTetrahedralCorotationalForceField<DataTypes>* ff;

};

topology::PointData<sofa::type::vector<Mat3x3> > pointInfo; ///< Internal point data
topology::EdgeData<sofa::type::vector<Mat3x3> > edgeInfo; ///< Internal edge data
topology::TetrahedronData<sofa::type::vector<TetrahedronRestInformation> > tetrahedronInfo; ///< Internal tetrahedron data

/** Method to initialize @sa TetrahedronRestInformation when a new Tetrahedron is created.
* Will be set as creation callback in the TetrahedronData @sa tetrahedronInfo
*/
void createTetrahedronRestInformation(Index, TetrahedronRestInformation& t,
const core::topology::BaseMeshTopology::Tetrahedron&,
const sofa::type::vector<Index>&,
const sofa::type::vector<double>&);

sofa::core::topology::BaseMeshTopology* m_topology;
VecCoord _initialPoints;///< the intial positions of the points
Expand Down Expand Up @@ -203,8 +185,6 @@ class FastTetrahedralCorotationalForceField : public core::behavior::ForceField<


protected :
FTCFTetrahedronHandler* tetrahedronHandler;

static void computeQRRotation( Mat3x3 &r, const Coord *dp);

topology::EdgeData<sofa::type::vector<Mat3x3> > &getEdgeInfo() {return edgeInfo;}
Expand Down
169 changes: 86 additions & 83 deletions modules/SofaMiscFem/src/SofaMiscFem/FastTetrahedralCorotationalForceField.inl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -36,104 +36,101 @@ namespace sofa::component::forcefield
using sofa::core::topology::edgesInTetrahedronArray;

template< class DataTypes>
void FastTetrahedralCorotationalForceField<DataTypes>::FTCFTetrahedronHandler::applyCreateFunction(Index tetrahedronIndex,
void FastTetrahedralCorotationalForceField<DataTypes>::createTetrahedronRestInformation(Index tetrahedronIndex,
TetrahedronRestInformation &my_tinfo,
const core::topology::BaseMeshTopology::Tetrahedron &,
const sofa::type::vector<Index> &,
const sofa::type::vector<double> &)
{
if (ff)
const std::vector< Tetrahedron > &tetrahedronArray=this->m_topology->getTetrahedra() ;
// const std::vector< Edge> &edgeArray=m_topology->getEdges() ;
unsigned int j,k,l,m,n;
typename DataTypes::Real lambda=getLambda();
typename DataTypes::Real mu=getMu();
typename DataTypes::Real volume,val;
typename DataTypes::Coord point[4]; //shapeVector[4];
const typename DataTypes::VecCoord restPosition=this->mstate->read(core::ConstVecCoordId::restPosition())->getValue();

///describe the indices of the 4 tetrahedron vertices
const Tetrahedron &t= tetrahedronArray[tetrahedronIndex];
// BaseMeshTopology::EdgesInTetrahedron te=m_topology->getEdgesInTetrahedron(tetrahedronIndex);


// store the point position
for(j=0; j<4; ++j)
point[j]=(restPosition)[t[j]];
/// compute 6 times the rest volume
volume=dot(cross(point[1]-point[0],point[2]-point[0]),point[0]-point[3]);
/// store the rest volume
my_tinfo.restVolume=volume/6;
mu*=fabs(volume)/6;
lambda*=fabs(volume)/6;

// store shape vectors at the rest configuration
for(j=0; j<4; ++j)
{
const std::vector< Tetrahedron > &tetrahedronArray=ff->m_topology->getTetrahedra() ;
// const std::vector< Edge> &edgeArray=ff->m_topology->getEdges() ;
unsigned int j,k,l,m,n;
typename DataTypes::Real lambda=ff->getLambda();
typename DataTypes::Real mu=ff->getMu();
typename DataTypes::Real volume,val;
typename DataTypes::Coord point[4]; //shapeVector[4];
const typename DataTypes::VecCoord restPosition=ff->mstate->read(core::ConstVecCoordId::restPosition())->getValue();

///describe the indices of the 4 tetrahedron vertices
const Tetrahedron &t= tetrahedronArray[tetrahedronIndex];
// BaseMeshTopology::EdgesInTetrahedron te=ff->m_topology->getEdgesInTetrahedron(tetrahedronIndex);


// store the point position
for(j=0; j<4; ++j)
point[j]=(restPosition)[t[j]];
/// compute 6 times the rest volume
volume=dot(cross(point[1]-point[0],point[2]-point[0]),point[0]-point[3]);
/// store the rest volume
my_tinfo.restVolume=volume/6;
mu*=fabs(volume)/6;
lambda*=fabs(volume)/6;

// store shape vectors at the rest configuration
for(j=0; j<4; ++j)
{
if ((j%2)==0)
my_tinfo.shapeVector[j]=cross(point[(j+2)%4] - point[(j+1)%4],point[(j+3)%4] - point[(j+1)%4])/volume;
else
my_tinfo.shapeVector[j]= -cross(point[(j+2)%4] - point[(j+1)%4],point[(j+3)%4] - point[(j+1)%4])/volume;
}
if ((j%2)==0)
my_tinfo.shapeVector[j]=cross(point[(j+2)%4] - point[(j+1)%4],point[(j+3)%4] - point[(j+1)%4])/volume;
else
my_tinfo.shapeVector[j]= -cross(point[(j+2)%4] - point[(j+1)%4],point[(j+3)%4] - point[(j+1)%4])/volume;
}

/// compute the vertex stiffness of the linear elastic material, needed for addKToMatrix
for(j=0; j<4; ++j)
/// compute the vertex stiffness of the linear elastic material, needed for addKToMatrix
for(j=0; j<4; ++j)
{
// the linear stiffness matrix using shape vectors and Lame coefficients
val=mu*dot(my_tinfo.shapeVector[j],my_tinfo.shapeVector[j]);
for(m=0; m<3; ++m)
{
// the linear stiffness matrix using shape vectors and Lame coefficients
val=mu*dot(my_tinfo.shapeVector[j],my_tinfo.shapeVector[j]);
for(m=0; m<3; ++m)
for(n=m; n<3; ++n)
{
for(n=m; n<3; ++n)
{
my_tinfo.linearDfDxDiag[j][m][n]=lambda*my_tinfo.shapeVector[j][n]*my_tinfo.shapeVector[j][m]+
mu*my_tinfo.shapeVector[j][n]*my_tinfo.shapeVector[j][m];
my_tinfo.linearDfDxDiag[j][m][n]=lambda*my_tinfo.shapeVector[j][n]*my_tinfo.shapeVector[j][m]+
mu*my_tinfo.shapeVector[j][n]*my_tinfo.shapeVector[j][m];

if (m==n)
{
my_tinfo.linearDfDxDiag[j][m][m]+=Real(val);
} else
my_tinfo.linearDfDxDiag[j][n][m]=my_tinfo.linearDfDxDiag[j][m][n];
}
if (m==n)
{
my_tinfo.linearDfDxDiag[j][m][m]+=Real(val);
} else
my_tinfo.linearDfDxDiag[j][n][m]=my_tinfo.linearDfDxDiag[j][m][n];
}
}
}

/// compute the edge stiffness of the linear elastic material
for(j=0; j<6; ++j)
{
core::topology::BaseMeshTopology::Edge e=ff->m_topology->getLocalEdgesInTetrahedron(j);
k=e[0];
l=e[1];
/// compute the edge stiffness of the linear elastic material
for(j=0; j<6; ++j)
{
core::topology::BaseMeshTopology::Edge e=this->m_topology->getLocalEdgesInTetrahedron(j);
k=e[0];
l=e[1];

// store the rest edge vector
my_tinfo.restEdgeVector[j]=point[l]-point[k];
// store the rest edge vector
my_tinfo.restEdgeVector[j]=point[l]-point[k];

// the linear stiffness matrix using shape vectors and Lame coefficients
val=mu*dot(my_tinfo.shapeVector[l],my_tinfo.shapeVector[k]);
for(m=0; m<3; ++m)
// the linear stiffness matrix using shape vectors and Lame coefficients
val=mu*dot(my_tinfo.shapeVector[l],my_tinfo.shapeVector[k]);
for(m=0; m<3; ++m)
{
for(n=0; n<3; ++n)
{
for(n=0; n<3; ++n)
{
my_tinfo.linearDfDx[j][m][n]=lambda*my_tinfo.shapeVector[k][n]*my_tinfo.shapeVector[l][m]+
mu*my_tinfo.shapeVector[l][n]*my_tinfo.shapeVector[k][m];
my_tinfo.linearDfDx[j][m][n]=lambda*my_tinfo.shapeVector[k][n]*my_tinfo.shapeVector[l][m]+
mu*my_tinfo.shapeVector[l][n]*my_tinfo.shapeVector[k][m];

if (m==n)
{
my_tinfo.linearDfDx[j][m][m]+=Real(val);
}
if (m==n)
{
my_tinfo.linearDfDx[j][m][m]+=Real(val);
}
}
}
if (ff->decompositionMethod==QR_DECOMPOSITION) {
// compute the rotation matrix of the initial tetrahedron for the QR decomposition
computeQRRotation(my_tinfo.restRotation,my_tinfo.restEdgeVector);
} else if (ff->decompositionMethod==POLAR_DECOMPOSITION_MODIFIED) {
Mat3x3 Transformation;
Transformation[0]=point[1]-point[0];
Transformation[1]=point[2]-point[0];
Transformation[2]=point[3]-point[0];
helper::Decompose<Real>::polarDecomposition( Transformation, my_tinfo.restRotation );
}
}
if (decompositionMethod==QR_DECOMPOSITION) {
// compute the rotation matrix of the initial tetrahedron for the QR decomposition
computeQRRotation(my_tinfo.restRotation,my_tinfo.restEdgeVector);
} else if (decompositionMethod==POLAR_DECOMPOSITION_MODIFIED) {
Mat3x3 Transformation;
Transformation[0]=point[1]-point[0];
Transformation[1]=point[2]-point[0];
Transformation[2]=point[3]-point[0];
helper::Decompose<Real>::polarDecomposition( Transformation, my_tinfo.restRotation );
}
}

Expand All @@ -155,14 +152,13 @@ template <class DataTypes> FastTetrahedralCorotationalForceField<DataTypes>::Fas
, drawColor3(initData(&drawColor3, sofa::type::RGBAColor(0.0f, 1.0f, 1.0f, 1.0f), "drawColor3", " draw color for faces 3"))
, drawColor4(initData(&drawColor4, sofa::type::RGBAColor(0.5f, 1.0f, 1.0f, 1.0f), "drawColor4", " draw color for faces 4"))
, l_topology(initLink("topology", "link to the topology container"))
, tetrahedronHandler(nullptr)
{
tetrahedronHandler = new FTCFTetrahedronHandler(this,&tetrahedronInfo);

}

template <class DataTypes> FastTetrahedralCorotationalForceField<DataTypes>::~FastTetrahedralCorotationalForceField()
{
if (tetrahedronHandler) delete tetrahedronHandler;

}

template <class DataTypes> void FastTetrahedralCorotationalForceField<DataTypes>::init()
Expand Down Expand Up @@ -234,12 +230,19 @@ template <class DataTypes> void FastTetrahedralCorotationalForceField<DataTypes>
/// initialize the data structure associated with each tetrahedron
for (Index i=0; i<m_topology->getNbTetrahedra(); ++i)
{
tetrahedronHandler->applyCreateFunction(i,tetrahedronInf[i],m_topology->getTetrahedron(i),
createTetrahedronRestInformation(i,tetrahedronInf[i],m_topology->getTetrahedron(i),
(const type::vector< Index > )0,
(const type::vector< double >)0);
}
/// set the call back function upon creation of a tetrahedron
tetrahedronInfo.createTopologyHandler(m_topology,tetrahedronHandler);
tetrahedronInfo.createTopologyHandler(m_topology);
tetrahedronInfo.setCreationCallback([this](Index tetrahedronIndex, TetrahedronRestInformation& tetraInfo,
const core::topology::BaseMeshTopology::Tetrahedron& tetra,
const sofa::type::vector< Index >& ancestors,
const sofa::type::vector< double >& coefs)
{
createTetrahedronRestInformation(tetrahedronIndex, tetraInfo, tetra, ancestors, coefs);
});
tetrahedronInfo.endEdit();

updateTopologyInfo=true;
Expand Down
21 changes: 7 additions & 14 deletions modules/SofaMiscFem/src/SofaMiscFem/QuadBendingFEMForceField.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -170,19 +170,13 @@ class QuadBendingFEMForceField : public core::behavior::ForceField<DataTypes>
topology::PointData<sofa::type::vector<VertexInformation> > vertexInfo; ///< Internal point data
topology::EdgeData<sofa::type::vector<EdgeInformation> > edgeInfo; ///< Internal edge data

class QuadHandler : public topology::TopologyDataHandler<core::topology::BaseMeshTopology::Quad,type::vector<QuadInformation> >
{
public:
QuadHandler(QuadBendingFEMForceField<DataTypes>* _ff, topology::QuadData<sofa::type::vector<QuadInformation> >* _data) :
topology::TopologyDataHandler<core::topology::BaseMeshTopology::Quad, sofa::type::vector<QuadInformation> >(_data), ff(_ff) {}

void applyCreateFunction(unsigned int quadIndex, QuadInformation& ,
const core::topology::BaseMeshTopology::Quad & t,
const sofa::type::vector< unsigned int > &,
const sofa::type::vector< double > &);
protected:
QuadBendingFEMForceField<DataTypes>* ff;
};
/** Method to initialize @sa QuadInformation when a new Quad is created.
* Will be set as creation callback in the QuadData @sa quadInfo
*/
void createQuadInformation(unsigned int quadIndex, QuadInformation&,
const core::topology::BaseMeshTopology::Quad& t,
const sofa::type::vector< unsigned int >&,
const sofa::type::vector< double >&);

sofa::core::topology::BaseMeshTopology* m_topology;

Expand Down Expand Up @@ -229,7 +223,6 @@ protected :
Data<type::vector<Real> > f_poisson; ///< Poisson ratio in Hooke's law (vector)
Data<type::vector<Real> > f_young; ///< Young modulus in Hooke's law (vector)
Data<Real> f_thickness;
QuadHandler* quadHandler;

/// Link to be set to the topology container in the component graph.
SingleLink<QuadBendingFEMForceField<DataTypes>, sofa::core::topology::BaseMeshTopology, BaseLink::FLAG_STOREPATH | BaseLink::FLAG_STRONGLINK> l_topology;
Expand Down
Loading