Skip to content

Commit

Permalink
[Mapping.NonLinear] A base class for (almost) all non linear mappings (…
Browse files Browse the repository at this point in the history 
…#5006)

* Introduce base class for non-linear mappings

* refactor DistanceFromTargetMapping

* refactor DistanceMapping

* refactor SquareDistanceMapping

* refactor SquareMapping

* delegate on updateK

* K can be private

* factorize baseMatrix

* baseMatrices can be private

* rename base class

* rename files

* fix

* second fix

* unused variables

* rename private members

* change protected member
  • Loading branch information
@alxbilger
alxbilger authored Oct 1, 2024
1 parent ede696a commit 465246f
Show file tree
Hide file tree
Showing 19 changed files with 651 additions and 959 deletions.
2 changes: 2 additions & 0 deletions Sofa/Component/Mapping/NonLinear/CMakeLists.txt
View file
Original file line number Diff line number Diff line change
Expand Up @@ -8,6 +8,8 @@ set(HEADER_FILES
${SOFACOMPONENTMAPPINGNONLINEAR_SOURCE_DIR}/init.h
${SOFACOMPONENTMAPPINGNONLINEAR_SOURCE_DIR}/AreaMapping.h
${SOFACOMPONENTMAPPINGNONLINEAR_SOURCE_DIR}/AreaMapping.inl
${SOFACOMPONENTMAPPINGNONLINEAR_SOURCE_DIR}/BaseNonLinearMapping.h
${SOFACOMPONENTMAPPINGNONLINEAR_SOURCE_DIR}/BaseNonLinearMapping.inl
${SOFACOMPONENTMAPPINGNONLINEAR_SOURCE_DIR}/DistanceFromTargetMapping.h
${SOFACOMPONENTMAPPINGNONLINEAR_SOURCE_DIR}/DistanceFromTargetMapping.inl
${SOFACOMPONENTMAPPINGNONLINEAR_SOURCE_DIR}/DistanceMapping.h
Expand Down
48 changes: 14 additions & 34 deletions Sofa/Component/Mapping/NonLinear/src/sofa/component/mapping/nonlinear/AreaMapping.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -22,7 +22,7 @@
#pragma once

#include <sofa/component/mapping/nonlinear/config.h>
#include <sofa/core/Mapping.h>
#include <sofa/component/mapping/nonlinear/BaseNonLinearMapping.h>
#include <sofa/component/mapping/nonlinear/NonLinearMappingData.h>
#include <sofa/core/topology/BaseMeshTopology.h>
#include <sofa/linearalgebra/EigenSparseMatrix.h>
Expand All @@ -32,17 +32,16 @@ namespace sofa::component::mapping::nonlinear
{

template <class TIn, class TOut>
class AreaMapping : public core::Mapping<TIn, TOut>, public StabilizedNonLinearMappingData
class AreaMapping : public BaseNonLinearMapping<TIn, TOut, true>
{
public:
SOFA_CLASS(SOFA_TEMPLATE2(AreaMapping,TIn,TOut), SOFA_TEMPLATE2(core::Mapping,TIn,TOut));
SOFA_CLASS(SOFA_TEMPLATE2(AreaMapping,TIn,TOut), SOFA_TEMPLATE3(BaseNonLinearMapping,TIn,TOut, true));

using In = TIn;
using Out = TOut;

using Real = Real_t<Out>;

typedef linearalgebra::EigenSparseMatrix<TIn,TOut> SparseMatrixEigen;
static constexpr auto Nin = In::deriv_total_size;

SingleLink<AreaMapping<TIn, TOut>, sofa::core::topology::BaseMeshTopology, BaseLink::FLAG_STOREPATH | BaseLink::FLAG_STRONGLINK> l_topology;
Expand All @@ -52,44 +51,25 @@ class AreaMapping : public core::Mapping<TIn, TOut>, public StabilizedNonLinearM
void init() override;

void apply(const core::MechanicalParams* mparams, DataVecCoord_t<Out>& out, const DataVecCoord_t<In>& in) override;
void applyJ(const core::MechanicalParams* mparams, DataVecDeriv_t<Out>& out, const DataVecDeriv_t<In>& in) override;
void applyJT(const core::MechanicalParams* mparams, DataVecDeriv_t<In>& out, const DataVecDeriv_t<Out>& in) override;
void applyJT(const core::ConstraintParams *cparams, DataMatrixDeriv_t<In>& out, const DataMatrixDeriv_t<Out>& in) override;
void applyDJT(const core::MechanicalParams* mparams, core::MultiVecDerivId parentForceId, core::ConstMultiVecDerivId childForceId) override;

void updateK( const core::MechanicalParams* mparams, core::ConstMultiVecDerivId childForceId) override;
const linearalgebra::BaseMatrix* getK() override;
void buildGeometricStiffnessMatrix(sofa::core::GeometricStiffnessMatrix* matrices) override;

const type::vector<sofa::linearalgebra::BaseMatrix*>* getJs() override;

protected:
AreaMapping();

using SparseKMatrixEigen = linearalgebra::EigenSparseMatrix<TIn,TIn>;

SparseMatrixEigen jacobian; ///< Jacobian of the mapping
type::vector<linearalgebra::BaseMatrix*> baseMatrices; ///< Jacobian of the mapping, in a vector
typename AreaMapping::SparseKMatrixEigen K; ///< Assembled geometric stiffness matrix

/**
* @brief Represents an entry in the Jacobian matrix.
*
* The JacobianEntry struct is used to store information about an entry in the
* Jacobian matrix, specifically the vertex identifier and the corresponding
* Jacobian value. It also provides a comparison operator for sorting entries
* by vertex ID.
*/
struct JacobianEntry
{
sofa::Index vertexId;
typename In::Coord jacobianValue;
bool operator<(const JacobianEntry& other) const { return vertexId < other.vertexId;}
};
void matrixFreeApplyDJT(const core::MechanicalParams* mparams, Real kFactor,
Data<VecDeriv_t<In> >& parentForce,
const Data<VecDeriv_t<In> >& parentDisplacement,
const Data<VecDeriv_t<Out> >& childForce) override;

const VecCoord_t<In>* m_vertices{nullptr};
using typename Inherit1::SparseKMatrixEigen;

void doUpdateK(
const core::MechanicalParams* mparams, const Data<VecDeriv_t<Out> >& childForce,
SparseKMatrixEigen& matrix) override;

const VecCoord_t<In>* m_vertices{nullptr};

using JacobianEntry = typename Inherit1::JacobianEntry;
};

#if !defined(SOFA_COMPONENT_MAPPING_NONLINEAR_AREAMAPPING_CPP)
Expand Down
158 changes: 45 additions & 113 deletions Sofa/Component/Mapping/NonLinear/src/sofa/component/mapping/nonlinear/AreaMapping.inl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -21,8 +21,8 @@
******************************************************************************/
#pragma once
#include <sofa/component/mapping/nonlinear/AreaMapping.h>
#include <sofa/component/mapping/nonlinear/BaseNonLinearMapping.inl>
#include <sofa/core/BaseLocalMappingMatrix.h>
#include <sofa/linearalgebra/CompressedRowSparseMatrixConstraintEigenUtils.h>

namespace sofa::component::mapping::nonlinear
{
Expand Down Expand Up @@ -123,10 +123,7 @@ void AreaMapping<TIn, TOut>::init()
typename core::behavior::MechanicalState<In>::ReadVecCoord pos = this->getFromModel()->readPositions();

this->getToModel()->resize( nbTriangles );
jacobian.resizeBlocks(nbTriangles, pos.size());

baseMatrices.resize( 1 );
baseMatrices[0] = &jacobian;
this->m_jacobian.resizeBlocks(nbTriangles, pos.size());

Inherit1::init();

Expand All @@ -140,14 +137,16 @@ template <class TIn, class TOut>
void AreaMapping<TIn, TOut>::apply(const core::MechanicalParams* mparams,
DataVecCoord_t<Out>& out, const DataVecCoord_t<In>& in)
{
SOFA_UNUSED( mparams );

helper::WriteOnlyAccessor< Data<VecCoord_t<Out>> > _out = out;
helper::ReadAccessor< Data<VecCoord_t<In> > > _in = in;

m_vertices = _in.operator->();

const auto& triangles = l_topology->getTriangles();

jacobian.clear();
this->m_jacobian.clear();

for (unsigned int triangleId = 0; triangleId < triangles.size(); ++triangleId)
{
Expand Down Expand Up @@ -177,132 +176,79 @@ void AreaMapping<TIn, TOut>::apply(const core::MechanicalParams* mparams,
//insertion in increasing column order
std::sort(jacobianEntries.begin(), jacobianEntries.end());

jacobian.beginRow(triangleId);
this->m_jacobian.beginRow(triangleId);
for (const auto& [vertexId, jacobianValue] : jacobianEntries)
{
for (unsigned d = 0; d < In::spatial_dimensions; ++d)
{
jacobian.insertBack(triangleId, vertexId * Nin + d, jacobianValue[d]);
this->m_jacobian.insertBack(triangleId, vertexId * Nin + d, jacobianValue[d]);
}
}
}

jacobian.compress();
}

template <class TIn, class TOut>
void AreaMapping<TIn, TOut>::applyJ(const core::MechanicalParams* mparams,
DataVecDeriv_t<Out>& out, const DataVecDeriv_t<In>& in)
{
if( jacobian.rowSize() )
{
auto dOutWa = sofa::helper::getWriteOnlyAccessor(out);
auto dInRa = sofa::helper::getReadAccessor(in);
jacobian.mult(dOutWa.wref(),dInRa.ref());
}
}

template <class TIn, class TOut>
void AreaMapping<TIn, TOut>::applyJT(const core::MechanicalParams* mparams,
DataVecDeriv_t<In>& out, const DataVecDeriv_t<Out>& in)
{
if( jacobian.rowSize() )
{
auto dOutRa = sofa::helper::getReadAccessor(in);
auto dInWa = sofa::helper::getWriteOnlyAccessor(out);
jacobian.addMultTranspose(dInWa.wref(),dOutRa.ref());
}
this->m_jacobian.compress();
}

template <class TIn, class TOut>
void AreaMapping<TIn, TOut>::applyJT(const core::ConstraintParams* cparams,
DataMatrixDeriv_t<In>& out, const DataMatrixDeriv_t<Out>& in)
void AreaMapping<TIn, TOut>::matrixFreeApplyDJT(
const core::MechanicalParams* mparams, Real kFactor,
Data<VecDeriv_t<In>>& parentForce,
const Data<VecDeriv_t<In>>& parentDisplacement,
const Data<VecDeriv_t<Out>>& childForce)
{
SOFA_UNUSED(cparams);
auto childMatRa = sofa::helper::getReadAccessor(in);
auto parentMatWa = sofa::helper::getWriteAccessor(out);
addMultTransposeEigen(parentMatWa.wref(), jacobian.compressedMatrix, childMatRa.ref());
}
SOFA_UNUSED(mparams);
const unsigned geometricStiffness = this->d_geometricStiffness.getValue().getSelectedId();

template <class TIn, class TOut>
void AreaMapping<TIn, TOut>::applyDJT(const core::MechanicalParams* mparams,
core::MultiVecDerivId parentForceId, core::ConstMultiVecDerivId childForceId)
{
if (!m_vertices)
{
return;
}
helper::WriteAccessor parentForceAccessor(parentForce);
helper::ReadAccessor parentDisplacementAccessor(parentDisplacement);
helper::ReadAccessor childForceAccessor(childForce);

const unsigned geometricStiffness = d_geometricStiffness.getValue().getSelectedId();
if( !geometricStiffness )
const auto& triangles = l_topology->getTriangles();
for (unsigned int triangleId = 0; triangleId < triangles.size(); ++triangleId)
{
return;
}

helper::WriteAccessor<Data<VecDeriv_t<In> > > parentForceAccessor(*parentForceId[this->fromModel.get()].write());
helper::ReadAccessor<Data<VecDeriv_t<In> > > parentDisplacementAccessor(*mparams->readDx(this->fromModel.get()));
const SReal kFactor = mparams->kFactor();
helper::ReadAccessor<Data<VecDeriv_t<Out> > > childForceAccessor(mparams->readF(this->toModel.get()));
const Deriv_t<Out>& childForceTri = childForceAccessor[triangleId];

if( K.compressedMatrix.nonZeros() )
{
K.addMult( parentForceAccessor.wref(), parentDisplacementAccessor.ref(), (typename In::Real)kFactor );
}
else
{
const auto& triangles = l_topology->getTriangles();
for (unsigned int triangleId = 0; triangleId < triangles.size(); ++triangleId)
if( childForceTri[0] < 0 || geometricStiffness==1 )
{
const Deriv_t<Out>& childForceTri = childForceAccessor[triangleId];

if( childForceTri[0] < 0 || geometricStiffness==1 )
{
const auto& triangle = triangles[triangleId];
const auto& triangle = triangles[triangleId];

const sofa::type::fixed_array<Coord_t<In>, 3> v{
(*m_vertices)[triangle[0]],
(*m_vertices)[triangle[1]],
(*m_vertices)[triangle[2]]
};
const sofa::type::fixed_array<Coord_t<In>, 3> v{
(*m_vertices)[triangle[0]],
(*m_vertices)[triangle[1]],
(*m_vertices)[triangle[2]]
};

//it's a 3x3 matrix, where each entry is a 3x3 matrix
const auto d2Area_d2x = computeSecondDerivativeArea(v);
//it's a 3x3 matrix, where each entry is a 3x3 matrix
const auto d2Area_d2x = computeSecondDerivativeArea(v);

for (unsigned int i = 0; i < 3; ++i)
for (unsigned int i = 0; i < 3; ++i)
{
for (unsigned int j = 0; j < 3; ++j)
{
for (unsigned int j = 0; j < 3; ++j)
{
parentForceAccessor[triangle[i]] +=
kFactor
* d2Area_d2x[i][j]
* parentDisplacementAccessor[triangle[j]]
* childForceTri[0];
}
parentForceAccessor[triangle[i]] +=
kFactor
* d2Area_d2x[i][j]
* parentDisplacementAccessor[triangle[j]]
* childForceTri[0];
}
}
}
}
}

template <class TIn, class TOut>
void AreaMapping<TIn, TOut>::updateK(const core::MechanicalParams* mparams,
core::ConstMultiVecDerivId childForceId)
void AreaMapping<TIn, TOut>::doUpdateK(const core::MechanicalParams* mparams,
const Data<VecDeriv_t<Out>>& childForce, SparseKMatrixEigen& matrix)
{
SOFA_UNUSED(mparams);
const unsigned geometricStiffness = d_geometricStiffness.getValue().getSelectedId();
if( !geometricStiffness ) { K.resize(0,0); return; }
const unsigned geometricStiffness = this->d_geometricStiffness.getValue().getSelectedId();

helper::ReadAccessor<Data<VecDeriv_t<Out> > > childForce( *childForceId[this->toModel.get()].read() );

{
unsigned int kSize = this->fromModel->getSize();
K.resizeBlocks(kSize, kSize);
}
const helper::ReadAccessor childForceAccessor(childForce);

const auto& triangles = l_topology->getTriangles();
for (unsigned int triangleId = 0; triangleId < triangles.size(); ++triangleId)
{
const Deriv_t<Out>& childForceTri = childForce[triangleId];
const Deriv_t<Out>& childForceTri = childForceAccessor[triangleId];

if( childForceTri[0] < 0 || geometricStiffness==1 )
{
Expand All @@ -321,26 +267,18 @@ void AreaMapping<TIn, TOut>::updateK(const core::MechanicalParams* mparams,
{
for (unsigned int j = 0; j < 3; ++j)
{
K.addBlock(triangle[i], triangle[j], d2Area_d2x[i][j] * childForceTri[0]);
matrix.addBlock(triangle[i], triangle[j], d2Area_d2x[i][j] * childForceTri[0]);
}
}
}
}

K.compress();
}

template <class TIn, class TOut>
const linearalgebra::BaseMatrix* AreaMapping<TIn, TOut>::getK()
{
return &K;
}

template <class TIn, class TOut>
void AreaMapping<TIn, TOut>::buildGeometricStiffnessMatrix(
sofa::core::GeometricStiffnessMatrix* matrices)
{
const unsigned& geometricStiffness = d_geometricStiffness.getValue().getSelectedId();
const unsigned& geometricStiffness = this->d_geometricStiffness.getValue().getSelectedId();
if( !geometricStiffness )
{
return;
Expand Down Expand Up @@ -379,11 +317,5 @@ void AreaMapping<TIn, TOut>::buildGeometricStiffnessMatrix(

}

template <class TIn, class TOut>
const type::vector<sofa::linearalgebra::BaseMatrix*>* AreaMapping<TIn, TOut>::
getJs()
{
return &baseMatrices;
}

}
Loading

0 comments on commit 465246f

Please sign in to comment.