From 6ba59228ef8eb221805d464ce2ed8e1b398bf47f Mon Sep 17 00:00:00 2001 From: Alex Bilger Date: Wed, 18 Sep 2024 14:49:12 +0200 Subject: [PATCH 01/16] Introduce base class for non-linear mappings --- .../Mapping/NonLinear/CMakeLists.txt | 3 + .../component/mapping/nonlinear/AreaMapping.h | 41 +---- .../mapping/nonlinear/AreaMapping.inl | 147 ++++++------------ .../nonlinear/AssembledNonLinearMapping.cpp | 33 ++++ .../nonlinear/AssembledNonLinearMapping.h | 91 +++++++++++ .../nonlinear/AssembledNonLinearMapping.inl | 118 ++++++++++++++ .../mapping/nonlinear/VolumeMapping.h | 41 +---- .../mapping/nonlinear/VolumeMapping.inl | 143 ++++++----------- ...essedRowSparseMatrixConstraintEigenUtils.h | 1 + 9 files changed, 351 insertions(+), 267 deletions(-) create mode 100644 Sofa/Component/Mapping/NonLinear/src/sofa/component/mapping/nonlinear/AssembledNonLinearMapping.cpp create mode 100644 Sofa/Component/Mapping/NonLinear/src/sofa/component/mapping/nonlinear/AssembledNonLinearMapping.h create mode 100644 Sofa/Component/Mapping/NonLinear/src/sofa/component/mapping/nonlinear/AssembledNonLinearMapping.inl diff --git a/Sofa/Component/Mapping/NonLinear/CMakeLists.txt b/Sofa/Component/Mapping/NonLinear/CMakeLists.txt index c5f5548c6ed..5a5ce42e3ab 100644 --- a/Sofa/Component/Mapping/NonLinear/CMakeLists.txt +++ b/Sofa/Component/Mapping/NonLinear/CMakeLists.txt @@ -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}/AssembledNonLinearMapping.h + ${SOFACOMPONENTMAPPINGNONLINEAR_SOURCE_DIR}/AssembledNonLinearMapping.inl ${SOFACOMPONENTMAPPINGNONLINEAR_SOURCE_DIR}/DistanceFromTargetMapping.h ${SOFACOMPONENTMAPPINGNONLINEAR_SOURCE_DIR}/DistanceFromTargetMapping.inl ${SOFACOMPONENTMAPPINGNONLINEAR_SOURCE_DIR}/DistanceMapping.h @@ -29,6 +31,7 @@ set(HEADER_FILES set(SOURCE_FILES ${SOFACOMPONENTMAPPINGNONLINEAR_SOURCE_DIR}/init.cpp ${SOFACOMPONENTMAPPINGNONLINEAR_SOURCE_DIR}/AreaMapping.cpp + ${SOFACOMPONENTMAPPINGNONLINEAR_SOURCE_DIR}/AssembledNonLinearMapping.cpp ${SOFACOMPONENTMAPPINGNONLINEAR_SOURCE_DIR}/DistanceFromTargetMapping.cpp ${SOFACOMPONENTMAPPINGNONLINEAR_SOURCE_DIR}/DistanceMapping.cpp ${SOFACOMPONENTMAPPINGNONLINEAR_SOURCE_DIR}/DistanceMultiMapping.cpp diff --git a/Sofa/Component/Mapping/NonLinear/src/sofa/component/mapping/nonlinear/AreaMapping.h b/Sofa/Component/Mapping/NonLinear/src/sofa/component/mapping/nonlinear/AreaMapping.h index 91f9cee3fc3..24457964f22 100644 --- a/Sofa/Component/Mapping/NonLinear/src/sofa/component/mapping/nonlinear/AreaMapping.h +++ b/Sofa/Component/Mapping/NonLinear/src/sofa/component/mapping/nonlinear/AreaMapping.h @@ -22,7 +22,7 @@ #pragma once #include -#include +#include #include #include #include @@ -32,17 +32,16 @@ namespace sofa::component::mapping::nonlinear { template -class AreaMapping : public core::Mapping, public StabilizedNonLinearMappingData +class AreaMapping : public AssembledNonLinearMapping { public: - SOFA_CLASS(SOFA_TEMPLATE2(AreaMapping,TIn,TOut), SOFA_TEMPLATE2(core::Mapping,TIn,TOut)); + SOFA_CLASS(SOFA_TEMPLATE2(AreaMapping,TIn,TOut), SOFA_TEMPLATE3(AssembledNonLinearMapping,TIn,TOut, true)); using In = TIn; using Out = TOut; using Real = Real_t; - typedef linearalgebra::EigenSparseMatrix SparseMatrixEigen; static constexpr auto Nin = In::deriv_total_size; SingleLink, sofa::core::topology::BaseMeshTopology, BaseLink::FLAG_STOREPATH | BaseLink::FLAG_STRONGLINK> l_topology; @@ -52,44 +51,20 @@ class AreaMapping : public core::Mapping, public StabilizedNonLinearM void init() override; void apply(const core::MechanicalParams* mparams, DataVecCoord_t& out, const DataVecCoord_t& in) override; - void applyJ(const core::MechanicalParams* mparams, DataVecDeriv_t& out, const DataVecDeriv_t& in) override; - void applyJT(const core::MechanicalParams* mparams, DataVecDeriv_t& out, const DataVecDeriv_t& in) override; - void applyJT(const core::ConstraintParams *cparams, DataMatrixDeriv_t& out, const DataMatrixDeriv_t& 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* getJs() override; - protected: AreaMapping(); - using SparseKMatrixEigen = linearalgebra::EigenSparseMatrix; - - SparseMatrixEigen jacobian; ///< Jacobian of the mapping - type::vector 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 >& parentForce, + const Data >& parentDisplacement, + const Data >& childForce) override; const VecCoord_t* m_vertices{nullptr}; - + using Inherit1::JacobianEntry; }; #if !defined(SOFA_COMPONENT_MAPPING_NONLINEAR_AREAMAPPING_CPP) diff --git a/Sofa/Component/Mapping/NonLinear/src/sofa/component/mapping/nonlinear/AreaMapping.inl b/Sofa/Component/Mapping/NonLinear/src/sofa/component/mapping/nonlinear/AreaMapping.inl index 88d62ca8b5d..bc3f2c741fe 100644 --- a/Sofa/Component/Mapping/NonLinear/src/sofa/component/mapping/nonlinear/AreaMapping.inl +++ b/Sofa/Component/Mapping/NonLinear/src/sofa/component/mapping/nonlinear/AreaMapping.inl @@ -21,8 +21,8 @@ ******************************************************************************/ #pragma once #include +#include #include -#include namespace sofa::component::mapping::nonlinear { @@ -123,10 +123,10 @@ void AreaMapping::init() typename core::behavior::MechanicalState::ReadVecCoord pos = this->getFromModel()->readPositions(); this->getToModel()->resize( nbTriangles ); - jacobian.resizeBlocks(nbTriangles, pos.size()); + this->jacobian.resizeBlocks(nbTriangles, pos.size()); - baseMatrices.resize( 1 ); - baseMatrices[0] = &jacobian; + this->baseMatrices.resize( 1 ); + this->baseMatrices[0] = &this->jacobian; Inherit1::init(); @@ -147,7 +147,7 @@ void AreaMapping::apply(const core::MechanicalParams* mparams, const auto& triangles = l_topology->getTriangles(); - jacobian.clear(); + this->jacobian.clear(); for (unsigned int triangleId = 0; triangleId < triangles.size(); ++triangleId) { @@ -177,107 +177,60 @@ void AreaMapping::apply(const core::MechanicalParams* mparams, //insertion in increasing column order std::sort(jacobianEntries.begin(), jacobianEntries.end()); - jacobian.beginRow(triangleId); + this->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->jacobian.insertBack(triangleId, vertexId * Nin + d, jacobianValue[d]); } } } - jacobian.compress(); + this->jacobian.compress(); } template -void AreaMapping::applyJ(const core::MechanicalParams* mparams, - DataVecDeriv_t& out, const DataVecDeriv_t& in) +void AreaMapping::matrixFreeApplyDJT( + const core::MechanicalParams* mparams, Real kFactor, + Data>& parentForce, + const Data>& parentDisplacement, + const Data>& childForce) { - if( jacobian.rowSize() ) - { - auto dOutWa = sofa::helper::getWriteOnlyAccessor(out); - auto dInRa = sofa::helper::getReadAccessor(in); - jacobian.mult(dOutWa.wref(),dInRa.ref()); - } -} - -template -void AreaMapping::applyJT(const core::MechanicalParams* mparams, - DataVecDeriv_t& out, const DataVecDeriv_t& in) -{ - if( jacobian.rowSize() ) - { - auto dOutRa = sofa::helper::getReadAccessor(in); - auto dInWa = sofa::helper::getWriteOnlyAccessor(out); - jacobian.addMultTranspose(dInWa.wref(),dOutRa.ref()); - } -} - -template -void AreaMapping::applyJT(const core::ConstraintParams* cparams, - DataMatrixDeriv_t& out, const DataMatrixDeriv_t& in) -{ - 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 -void AreaMapping::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 > > parentForceAccessor(*parentForceId[this->fromModel.get()].write()); - helper::ReadAccessor > > parentDisplacementAccessor(*mparams->readDx(this->fromModel.get())); - const SReal kFactor = mparams->kFactor(); - helper::ReadAccessor > > childForceAccessor(mparams->readF(this->toModel.get())); + const Deriv_t& 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& childForceTri = childForceAccessor[triangleId]; - - if( childForceTri[0] < 0 || geometricStiffness==1 ) - { - const auto& triangle = triangles[triangleId]; + const auto& triangle = triangles[triangleId]; - const sofa::type::fixed_array, 3> v{ - (*m_vertices)[triangle[0]], - (*m_vertices)[triangle[1]], - (*m_vertices)[triangle[2]] - }; + const sofa::type::fixed_array, 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]; } } } @@ -289,14 +242,14 @@ void AreaMapping::updateK(const core::MechanicalParams* mparams, core::ConstMultiVecDerivId childForceId) { 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(); + if( !geometricStiffness ) { this->K.resize(0,0); return; } helper::ReadAccessor > > childForce( *childForceId[this->toModel.get()].read() ); { unsigned int kSize = this->fromModel->getSize(); - K.resizeBlocks(kSize, kSize); + this->K.resizeBlocks(kSize, kSize); } const auto& triangles = l_topology->getTriangles(); @@ -321,26 +274,20 @@ void AreaMapping::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]); + this->K.addBlock(triangle[i], triangle[j], d2Area_d2x[i][j] * childForceTri[0]); } } } } - K.compress(); -} - -template -const linearalgebra::BaseMatrix* AreaMapping::getK() -{ - return &K; + this->K.compress(); } template void AreaMapping::buildGeometricStiffnessMatrix( sofa::core::GeometricStiffnessMatrix* matrices) { - const unsigned& geometricStiffness = d_geometricStiffness.getValue().getSelectedId(); + const unsigned& geometricStiffness = this->d_geometricStiffness.getValue().getSelectedId(); if( !geometricStiffness ) { return; @@ -379,11 +326,5 @@ void AreaMapping::buildGeometricStiffnessMatrix( } -template -const type::vector* AreaMapping:: -getJs() -{ - return &baseMatrices; -} } diff --git a/Sofa/Component/Mapping/NonLinear/src/sofa/component/mapping/nonlinear/AssembledNonLinearMapping.cpp b/Sofa/Component/Mapping/NonLinear/src/sofa/component/mapping/nonlinear/AssembledNonLinearMapping.cpp new file mode 100644 index 00000000000..fcc81791fc5 --- /dev/null +++ b/Sofa/Component/Mapping/NonLinear/src/sofa/component/mapping/nonlinear/AssembledNonLinearMapping.cpp @@ -0,0 +1,33 @@ +/****************************************************************************** +* SOFA, Simulation Open-Framework Architecture * +* (c) 2006 INRIA, USTL, UJF, CNRS, MGH * +* * +* This program is free software; you can redistribute it and/or modify it * +* under the terms of the GNU Lesser General Public License as published by * +* the Free Software Foundation; either version 2.1 of the License, or (at * +* your option) any later version. * +* * +* This program is distributed in the hope that it will be useful, but WITHOUT * +* ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or * +* FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License * +* for more details. * +* * +* You should have received a copy of the GNU Lesser General Public License * +* along with this program. If not, see . * +******************************************************************************* +* Authors: The SOFA Team and external contributors (see Authors.txt) * +* * +* Contact information: contact@sofa-framework.org * +******************************************************************************/ +#define SOFA_COMPONENT_MAPPING_NONLINEAR_ASSEMBLEDNONLINEARMAPPING_CPP +#include +#include + +namespace sofa::component::mapping::nonlinear +{ + +using namespace defaulttype; + +template class SOFA_COMPONENT_MAPPING_NONLINEAR_API AssembledNonLinearMapping< Vec3Types, Vec1Types, true >; + +} diff --git a/Sofa/Component/Mapping/NonLinear/src/sofa/component/mapping/nonlinear/AssembledNonLinearMapping.h b/Sofa/Component/Mapping/NonLinear/src/sofa/component/mapping/nonlinear/AssembledNonLinearMapping.h new file mode 100644 index 00000000000..d91ecf1ff56 --- /dev/null +++ b/Sofa/Component/Mapping/NonLinear/src/sofa/component/mapping/nonlinear/AssembledNonLinearMapping.h @@ -0,0 +1,91 @@ +/****************************************************************************** +* SOFA, Simulation Open-Framework Architecture * +* (c) 2006 INRIA, USTL, UJF, CNRS, MGH * +* * +* This program is free software; you can redistribute it and/or modify it * +* under the terms of the GNU Lesser General Public License as published by * +* the Free Software Foundation; either version 2.1 of the License, or (at * +* your option) any later version. * +* * +* This program is distributed in the hope that it will be useful, but WITHOUT * +* ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or * +* FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License * +* for more details. * +* * +* You should have received a copy of the GNU Lesser General Public License * +* along with this program. If not, see . * +******************************************************************************* +* Authors: The SOFA Team and external contributors (see Authors.txt) * +* * +* Contact information: contact@sofa-framework.org * +******************************************************************************/ +#pragma once +#include +#include +#include +#include + + +namespace sofa::component::mapping::nonlinear +{ + +template +class AssembledNonLinearMapping : public core::Mapping, public NonLinearMappingData +{ +public: + SOFA_CLASS( + SOFA_TEMPLATE3(AssembledNonLinearMapping,TIn,TOut, HasStabilizedGeometricStiffness), + SOFA_TEMPLATE2(core::Mapping,TIn,TOut)); + + using In = TIn; + using Out = TOut; + + using Real = Real_t; + + typedef linearalgebra::EigenSparseMatrix SparseMatrixEigen; + static constexpr auto Nin = In::deriv_total_size; + + void init() override; + + void applyJ(const core::MechanicalParams* mparams, DataVecDeriv_t& out, const DataVecDeriv_t& in) final; + void applyJT(const core::MechanicalParams* mparams, DataVecDeriv_t& out, const DataVecDeriv_t& in) final; + void applyJT(const core::ConstraintParams *cparams, DataMatrixDeriv_t& out, const DataMatrixDeriv_t& in) final; + void applyDJT(const core::MechanicalParams* mparams, core::MultiVecDerivId parentForceId, core::ConstMultiVecDerivId childForceId) final; + + const linearalgebra::BaseMatrix* getK() final; + const type::vector* getJs() override; + +protected: + virtual void matrixFreeApplyDJT(const core::MechanicalParams* mparams, + Real kFactor, + Data >& parentForce, + const Data >& parentDisplacement, + const Data >& childForce) = 0; + + using SparseKMatrixEigen = linearalgebra::EigenSparseMatrix; + + SparseMatrixEigen jacobian; ///< Jacobian of the mapping + type::vector baseMatrices; ///< Jacobian of the mapping, in a vector + 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;} + }; +}; + +#if !defined(SOFA_COMPONENT_MAPPING_NONLINEAR_ASSEMBLEDNONLINEARMAPPING_CPP) +extern template class SOFA_COMPONENT_MAPPING_NONLINEAR_API AssembledNonLinearMapping; +#endif + +} diff --git a/Sofa/Component/Mapping/NonLinear/src/sofa/component/mapping/nonlinear/AssembledNonLinearMapping.inl b/Sofa/Component/Mapping/NonLinear/src/sofa/component/mapping/nonlinear/AssembledNonLinearMapping.inl new file mode 100644 index 00000000000..d636d1ccad7 --- /dev/null +++ b/Sofa/Component/Mapping/NonLinear/src/sofa/component/mapping/nonlinear/AssembledNonLinearMapping.inl @@ -0,0 +1,118 @@ +/****************************************************************************** +* SOFA, Simulation Open-Framework Architecture * +* (c) 2006 INRIA, USTL, UJF, CNRS, MGH * +* * +* This program is free software; you can redistribute it and/or modify it * +* under the terms of the GNU Lesser General Public License as published by * +* the Free Software Foundation; either version 2.1 of the License, or (at * +* your option) any later version. * +* * +* This program is distributed in the hope that it will be useful, but WITHOUT * +* ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or * +* FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License * +* for more details. * +* * +* You should have received a copy of the GNU Lesser General Public License * +* along with this program. If not, see . * +******************************************************************************* +* Authors: The SOFA Team and external contributors (see Authors.txt) * +* * +* Contact information: contact@sofa-framework.org * +******************************************************************************/ +#pragma once + +#include +#include + +namespace sofa::component::mapping::nonlinear +{ + +template +void AssembledNonLinearMapping:: +init() +{ + core::Mapping::init(); +} + +template +void AssembledNonLinearMapping::applyJ( + const core::MechanicalParams* mparams, DataVecDeriv_t& out, + const DataVecDeriv_t& in) +{ + if( jacobian.rowSize() ) + { + auto dOutWa = sofa::helper::getWriteOnlyAccessor(out); + auto dInRa = sofa::helper::getReadAccessor(in); + jacobian.mult(dOutWa.wref(),dInRa.ref()); + } +} + +template +void AssembledNonLinearMapping::applyJT( + const core::MechanicalParams* mparams, DataVecDeriv_t& out, + const DataVecDeriv_t& in) +{ + if( jacobian.rowSize() ) + { + auto dOutRa = sofa::helper::getReadAccessor(in); + auto dInWa = sofa::helper::getWriteOnlyAccessor(out); + jacobian.addMultTranspose(dInWa.wref(),dOutRa.ref()); + } +} + +template +void AssembledNonLinearMapping::applyJT( + const core::ConstraintParams* cparams, DataMatrixDeriv_t& out, + const DataMatrixDeriv_t& in) +{ + SOFA_UNUSED(cparams); + auto childMatRa = sofa::helper::getReadAccessor(in); + auto parentMatWa = sofa::helper::getWriteAccessor(out); + addMultTransposeEigen(parentMatWa.wref(), jacobian.compressedMatrix, childMatRa.ref()); +} + +template +void AssembledNonLinearMapping::applyDJT( + const core::MechanicalParams* mparams, core::MultiVecDerivId parentForceId, + core::ConstMultiVecDerivId childForceId) +{ + const unsigned geometricStiffness = this->d_geometricStiffness.getValue().getSelectedId(); + if( !geometricStiffness ) + { + return; + } + + Data >& parentForce = *parentForceId[this->fromModel.get()].write(); + const Data >& parentDisplacement = *mparams->readDx(this->fromModel.get()); + + const SReal kFactor = mparams->kFactor(); + + if( K.compressedMatrix.nonZeros() ) + { + helper::WriteAccessor parentForceAccessor(parentForce); + helper::ReadAccessor parentDisplacementAccessor(parentDisplacement); + K.addMult( parentForceAccessor.wref(), parentDisplacementAccessor.ref(), static_cast(kFactor) ); + } + else + { + const Data >& childForce = *mparams->readF(this->toModel.get()); + + matrixFreeApplyDJT(mparams, static_cast(kFactor), + parentForce, parentDisplacement, childForce); + } +} + +template +const linearalgebra::BaseMatrix* AssembledNonLinearMapping::getK() +{ + return &K; +} + +template +const type::vector* AssembledNonLinearMapping::getJs() +{ + return &baseMatrices; +} + + +} diff --git a/Sofa/Component/Mapping/NonLinear/src/sofa/component/mapping/nonlinear/VolumeMapping.h b/Sofa/Component/Mapping/NonLinear/src/sofa/component/mapping/nonlinear/VolumeMapping.h index 80aae2b96d8..e677bb70503 100644 --- a/Sofa/Component/Mapping/NonLinear/src/sofa/component/mapping/nonlinear/VolumeMapping.h +++ b/Sofa/Component/Mapping/NonLinear/src/sofa/component/mapping/nonlinear/VolumeMapping.h @@ -22,7 +22,7 @@ #pragma once #include -#include +#include #include #include #include @@ -32,17 +32,16 @@ namespace sofa::component::mapping::nonlinear { template -class VolumeMapping : public core::Mapping, public StabilizedNonLinearMappingData +class VolumeMapping : public AssembledNonLinearMapping { public: - SOFA_CLASS(SOFA_TEMPLATE2(VolumeMapping,TIn,TOut), SOFA_TEMPLATE2(core::Mapping,TIn,TOut)); + SOFA_CLASS(SOFA_TEMPLATE2(VolumeMapping,TIn,TOut), SOFA_TEMPLATE3(AssembledNonLinearMapping,TIn,TOut, true)); using In = TIn; using Out = TOut; using Real = Real_t; - typedef linearalgebra::EigenSparseMatrix SparseMatrixEigen; static constexpr auto Nin = In::deriv_total_size; SingleLink, sofa::core::topology::BaseMeshTopology, BaseLink::FLAG_STOREPATH | BaseLink::FLAG_STRONGLINK> l_topology; @@ -52,44 +51,20 @@ class VolumeMapping : public core::Mapping, public StabilizedNonLinea void init() override; void apply(const core::MechanicalParams* mparams, DataVecCoord_t& out, const DataVecCoord_t& in) override; - void applyJ(const core::MechanicalParams* mparams, DataVecDeriv_t& out, const DataVecDeriv_t& in) override; - void applyJT(const core::MechanicalParams* mparams, DataVecDeriv_t& out, const DataVecDeriv_t& in) override; - void applyJT(const core::ConstraintParams *cparams, DataMatrixDeriv_t& out, const DataMatrixDeriv_t& 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* getJs() override; - protected: VolumeMapping(); - using SparseKMatrixEigen = linearalgebra::EigenSparseMatrix; - - SparseMatrixEigen jacobian; ///< Jacobian of the mapping - type::vector baseMatrices; ///< Jacobian of the mapping, in a vector - typename VolumeMapping::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 >& parentForce, + const Data >& parentDisplacement, + const Data >& childForce) override; const VecCoord_t* m_vertices{nullptr}; - + using Inherit1::JacobianEntry; }; #if !defined(SOFA_COMPONENT_MAPPING_NONLINEAR_VOLUMEMAPPING_CPP) diff --git a/Sofa/Component/Mapping/NonLinear/src/sofa/component/mapping/nonlinear/VolumeMapping.inl b/Sofa/Component/Mapping/NonLinear/src/sofa/component/mapping/nonlinear/VolumeMapping.inl index 4e7fea16004..94e6d4b266f 100644 --- a/Sofa/Component/Mapping/NonLinear/src/sofa/component/mapping/nonlinear/VolumeMapping.inl +++ b/Sofa/Component/Mapping/NonLinear/src/sofa/component/mapping/nonlinear/VolumeMapping.inl @@ -21,8 +21,8 @@ ******************************************************************************/ #pragma once #include +#include #include -#include namespace sofa::component::mapping::nonlinear { @@ -102,10 +102,10 @@ void VolumeMapping::init() typename core::behavior::MechanicalState::ReadVecCoord pos = this->getFromModel()->readPositions(); this->getToModel()->resize( nbTetrahedra ); - jacobian.resizeBlocks(nbTetrahedra, pos.size()); + this->jacobian.resizeBlocks(nbTetrahedra, pos.size()); - baseMatrices.resize( 1 ); - baseMatrices[0] = &jacobian; + this->baseMatrices.resize( 1 ); + this->baseMatrices[0] = &this->jacobian; Inherit1::init(); @@ -126,7 +126,7 @@ void VolumeMapping::apply(const core::MechanicalParams* mparams, const auto& tetrahedra = l_topology->getTetrahedra(); - jacobian.clear(); + this->jacobian.clear(); for (unsigned int tetId = 0; tetId < tetrahedra.size(); ++tetId) { @@ -166,115 +166,74 @@ void VolumeMapping::apply(const core::MechanicalParams* mparams, //insertion in increasing column order std::sort(jacobianEntries.begin(), jacobianEntries.end()); - jacobian.beginRow(tetId); + this->jacobian.beginRow(tetId); for (const auto& [vertexId, jacobianValue] : jacobianEntries) { for (unsigned d = 0; d < In::spatial_dimensions; ++d) { - jacobian.insertBack(tetId, vertexId * Nin + d, jacobianValue[d]); + this->jacobian.insertBack(tetId, vertexId * Nin + d, jacobianValue[d]); } } } - jacobian.compress(); + this->jacobian.compress(); } template -void VolumeMapping::applyJ(const core::MechanicalParams* mparams, - DataVecDeriv_t& out, const DataVecDeriv_t& in) +void VolumeMapping::matrixFreeApplyDJT( + const core::MechanicalParams* mparams, Real kFactor, + Data>& parentForce, + const Data>& parentDisplacement, + const Data>& childForce) { - if( jacobian.rowSize() ) - { - auto dOutWa = sofa::helper::getWriteOnlyAccessor(out); - auto dInRa = sofa::helper::getReadAccessor(in); - jacobian.mult(dOutWa.wref(),dInRa.ref()); - } -} - -template -void VolumeMapping::applyJT(const core::MechanicalParams* mparams, - DataVecDeriv_t& out, const DataVecDeriv_t& in) -{ - if( jacobian.rowSize() ) - { - auto dOutRa = sofa::helper::getReadAccessor(in); - auto dInWa = sofa::helper::getWriteOnlyAccessor(out); - jacobian.addMultTranspose(dInWa.wref(),dOutRa.ref()); - } -} - -template -void VolumeMapping::applyJT(const core::ConstraintParams* cparams, - DataMatrixDeriv_t& out, const DataMatrixDeriv_t& in) -{ - 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); -template -void VolumeMapping::applyDJT(const core::MechanicalParams* mparams, - core::MultiVecDerivId parentForceId, core::ConstMultiVecDerivId childForceId) -{ if (!m_vertices) { return; } - const unsigned geometricStiffness = d_geometricStiffness.getValue().getSelectedId(); - if( !geometricStiffness ) - { - return; - } + const unsigned geometricStiffness = this->d_geometricStiffness.getValue().getSelectedId(); - helper::WriteAccessor > > parentForceAccessor(*parentForceId[this->fromModel.get()].write()); - helper::ReadAccessor > > parentDisplacementAccessor(*mparams->readDx(this->fromModel.get())); - const SReal kFactor = mparams->kFactor(); - helper::ReadAccessor > > childForceAccessor(mparams->readF(this->toModel.get())); + helper::WriteAccessor parentForceAccessor(parentForce); + helper::ReadAccessor parentDisplacementAccessor(parentDisplacement); + helper::ReadAccessor childForceAccessor(childForce); - if( K.compressedMatrix.nonZeros() ) - { - K.addMult( parentForceAccessor.wref(), parentDisplacementAccessor.ref(), (typename In::Real)kFactor ); - } - else + const auto& tetrahedra = l_topology->getTetrahedra(); + for (unsigned int tetId = 0; tetId < tetrahedra.size(); ++tetId) { - const auto& tetrahedra = l_topology->getTetrahedra(); - for (unsigned int tetId = 0; tetId < tetrahedra.size(); ++tetId) - { - const Deriv_t& childForceTetra = childForceAccessor[tetId]; + const Deriv_t& childForceTetra = childForceAccessor[tetId]; - if( childForceTetra[0] < 0 || geometricStiffness==1 ) - { - const auto& tetra = tetrahedra[tetId]; + if( childForceTetra[0] < 0 || geometricStiffness==1 ) + { + const auto& tetra = tetrahedra[tetId]; - const sofa::type::fixed_array, 4> v{ - (*m_vertices)[tetra[0]], - (*m_vertices)[tetra[1]], - (*m_vertices)[tetra[2]], - (*m_vertices)[tetra[3]], - }; + const sofa::type::fixed_array, 4> v{ + (*m_vertices)[tetra[0]], + (*m_vertices)[tetra[1]], + (*m_vertices)[tetra[2]], + (*m_vertices)[tetra[3]], + }; - //it's a 4x4 matrix, where each entry is a 3x3 matrix - const auto d2Vol_d2x = computeSecondDerivativeVolume(v); + //it's a 4x4 matrix, where each entry is a 3x3 matrix + const auto d2Vol_d2x = computeSecondDerivativeVolume(v); - for (unsigned int i = 0; i < 4; ++i) + for (unsigned int i = 0; i < 4; ++i) + { + for (unsigned int j = i + 1; j < 4; ++j) //diagonal terms are omitted because they are null { - for (unsigned int j = i + 1; j < 4; ++j) //diagonal terms are omitted because they are null - { - parentForceAccessor[tetra[i]] += + parentForceAccessor[tetra[i]] += kFactor * d2Vol_d2x(i, j) * parentDisplacementAccessor[tetra[j]] * childForceTetra[0]; - //transpose - parentForceAccessor[tetra[j]] += + //transpose + parentForceAccessor[tetra[j]] += kFactor * d2Vol_d2x(j, i) * parentDisplacementAccessor[tetra[i]] * childForceTetra[0]; - } } } } @@ -286,14 +245,14 @@ void VolumeMapping::updateK(const core::MechanicalParams* mparams, core::ConstMultiVecDerivId childForceId) { 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(); + if( !geometricStiffness ) { this->K.resize(0,0); return; } helper::ReadAccessor > > childForce( *childForceId[this->toModel.get()].read() ); { unsigned int kSize = this->fromModel->getSize(); - K.resizeBlocks(kSize, kSize); + this->K.resizeBlocks(kSize, kSize); } const auto& tetrahedra = l_topology->getTetrahedra(); @@ -319,27 +278,21 @@ void VolumeMapping::updateK(const core::MechanicalParams* mparams, { for (unsigned int j = i+1; j < 4; ++j) //diagonal terms are omitted because they are null { - K.addBlock(tetra[i], tetra[j], d2Volume_d2x(i, j) * childForceTri[0]); - K.addBlock(tetra[j], tetra[i], d2Volume_d2x(j, i) * childForceTri[0]); + this->K.addBlock(tetra[i], tetra[j], d2Volume_d2x(i, j) * childForceTri[0]); + this->K.addBlock(tetra[j], tetra[i], d2Volume_d2x(j, i) * childForceTri[0]); } } } } - K.compress(); -} - -template -const linearalgebra::BaseMatrix* VolumeMapping::getK() -{ - return &K; + this->K.compress(); } template void VolumeMapping::buildGeometricStiffnessMatrix( sofa::core::GeometricStiffnessMatrix* matrices) { - const unsigned& geometricStiffness = d_geometricStiffness.getValue().getSelectedId(); + const unsigned& geometricStiffness = this->d_geometricStiffness.getValue().getSelectedId(); if( !geometricStiffness ) { return; @@ -380,11 +333,5 @@ void VolumeMapping::buildGeometricStiffnessMatrix( } -template -const type::vector* VolumeMapping:: -getJs() -{ - return &baseMatrices; -} } diff --git a/Sofa/framework/LinearAlgebra/src/sofa/linearalgebra/CompressedRowSparseMatrixConstraintEigenUtils.h b/Sofa/framework/LinearAlgebra/src/sofa/linearalgebra/CompressedRowSparseMatrixConstraintEigenUtils.h index dfdb2afd0c8..88da4face8e 100644 --- a/Sofa/framework/LinearAlgebra/src/sofa/linearalgebra/CompressedRowSparseMatrixConstraintEigenUtils.h +++ b/Sofa/framework/LinearAlgebra/src/sofa/linearalgebra/CompressedRowSparseMatrixConstraintEigenUtils.h @@ -19,6 +19,7 @@ * * * Contact information: contact@sofa-framework.org * ******************************************************************************/ +#pragma once #include #include #include From 8e956547a30afa1cced0ad3c2649dfd1a56a688d Mon Sep 17 00:00:00 2001 From: Alex Bilger Date: Wed, 18 Sep 2024 15:07:02 +0200 Subject: [PATCH 02/16] refactor DistanceFromTargetMapping --- .../Mapping/NonLinear/CMakeLists.txt | 1 - .../nonlinear/AssembledNonLinearMapping.cpp | 33 ----- .../nonlinear/AssembledNonLinearMapping.h | 4 - .../nonlinear/DistanceFromTargetMapping.h | 58 +++------ .../nonlinear/DistanceFromTargetMapping.inl | 123 +++++------------- 5 files changed, 54 insertions(+), 165 deletions(-) delete mode 100644 Sofa/Component/Mapping/NonLinear/src/sofa/component/mapping/nonlinear/AssembledNonLinearMapping.cpp diff --git a/Sofa/Component/Mapping/NonLinear/CMakeLists.txt b/Sofa/Component/Mapping/NonLinear/CMakeLists.txt index 5a5ce42e3ab..03545b647ce 100644 --- a/Sofa/Component/Mapping/NonLinear/CMakeLists.txt +++ b/Sofa/Component/Mapping/NonLinear/CMakeLists.txt @@ -31,7 +31,6 @@ set(HEADER_FILES set(SOURCE_FILES ${SOFACOMPONENTMAPPINGNONLINEAR_SOURCE_DIR}/init.cpp ${SOFACOMPONENTMAPPINGNONLINEAR_SOURCE_DIR}/AreaMapping.cpp - ${SOFACOMPONENTMAPPINGNONLINEAR_SOURCE_DIR}/AssembledNonLinearMapping.cpp ${SOFACOMPONENTMAPPINGNONLINEAR_SOURCE_DIR}/DistanceFromTargetMapping.cpp ${SOFACOMPONENTMAPPINGNONLINEAR_SOURCE_DIR}/DistanceMapping.cpp ${SOFACOMPONENTMAPPINGNONLINEAR_SOURCE_DIR}/DistanceMultiMapping.cpp diff --git a/Sofa/Component/Mapping/NonLinear/src/sofa/component/mapping/nonlinear/AssembledNonLinearMapping.cpp b/Sofa/Component/Mapping/NonLinear/src/sofa/component/mapping/nonlinear/AssembledNonLinearMapping.cpp deleted file mode 100644 index fcc81791fc5..00000000000 --- a/Sofa/Component/Mapping/NonLinear/src/sofa/component/mapping/nonlinear/AssembledNonLinearMapping.cpp +++ /dev/null @@ -1,33 +0,0 @@ -/****************************************************************************** -* SOFA, Simulation Open-Framework Architecture * -* (c) 2006 INRIA, USTL, UJF, CNRS, MGH * -* * -* This program is free software; you can redistribute it and/or modify it * -* under the terms of the GNU Lesser General Public License as published by * -* the Free Software Foundation; either version 2.1 of the License, or (at * -* your option) any later version. * -* * -* This program is distributed in the hope that it will be useful, but WITHOUT * -* ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or * -* FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License * -* for more details. * -* * -* You should have received a copy of the GNU Lesser General Public License * -* along with this program. If not, see . * -******************************************************************************* -* Authors: The SOFA Team and external contributors (see Authors.txt) * -* * -* Contact information: contact@sofa-framework.org * -******************************************************************************/ -#define SOFA_COMPONENT_MAPPING_NONLINEAR_ASSEMBLEDNONLINEARMAPPING_CPP -#include -#include - -namespace sofa::component::mapping::nonlinear -{ - -using namespace defaulttype; - -template class SOFA_COMPONENT_MAPPING_NONLINEAR_API AssembledNonLinearMapping< Vec3Types, Vec1Types, true >; - -} diff --git a/Sofa/Component/Mapping/NonLinear/src/sofa/component/mapping/nonlinear/AssembledNonLinearMapping.h b/Sofa/Component/Mapping/NonLinear/src/sofa/component/mapping/nonlinear/AssembledNonLinearMapping.h index d91ecf1ff56..3c7845fa129 100644 --- a/Sofa/Component/Mapping/NonLinear/src/sofa/component/mapping/nonlinear/AssembledNonLinearMapping.h +++ b/Sofa/Component/Mapping/NonLinear/src/sofa/component/mapping/nonlinear/AssembledNonLinearMapping.h @@ -84,8 +84,4 @@ class AssembledNonLinearMapping : public core::Mapping, public NonLin }; }; -#if !defined(SOFA_COMPONENT_MAPPING_NONLINEAR_ASSEMBLEDNONLINEARMAPPING_CPP) -extern template class SOFA_COMPONENT_MAPPING_NONLINEAR_API AssembledNonLinearMapping; -#endif - } diff --git a/Sofa/Component/Mapping/NonLinear/src/sofa/component/mapping/nonlinear/DistanceFromTargetMapping.h b/Sofa/Component/Mapping/NonLinear/src/sofa/component/mapping/nonlinear/DistanceFromTargetMapping.h index 17942e099c5..2218a80f9f2 100644 --- a/Sofa/Component/Mapping/NonLinear/src/sofa/component/mapping/nonlinear/DistanceFromTargetMapping.h +++ b/Sofa/Component/Mapping/NonLinear/src/sofa/component/mapping/nonlinear/DistanceFromTargetMapping.h @@ -22,8 +22,7 @@ #pragma once #include - -#include +#include #include #include #include @@ -61,28 +60,20 @@ struct BaseDistanceFromTargetMapping * @author Francois Faure */ template -class DistanceFromTargetMapping : public core::Mapping, public BaseDistanceFromTargetMapping, public NonLinearMappingData +class DistanceFromTargetMapping : public AssembledNonLinearMapping, public BaseDistanceFromTargetMapping { public: - SOFA_CLASS(SOFA_TEMPLATE2(DistanceFromTargetMapping,TIn,TOut), SOFA_TEMPLATE2(core::Mapping,TIn,TOut)); - - typedef core::Mapping Inherit; - typedef TIn In; - typedef TOut Out; - typedef typename Out::VecCoord OutVecCoord; - typedef typename Out::VecDeriv OutVecDeriv; - typedef typename Out::Coord OutCoord; - typedef typename Out::Deriv OutDeriv; - typedef typename Out::MatrixDeriv OutMatrixDeriv; - typedef typename Out::Real Real; - typedef typename In::Deriv InDeriv; - typedef typename In::MatrixDeriv InMatrixDeriv; - typedef typename In::Coord InCoord; - typedef typename In::VecCoord InVecCoord; - typedef typename In::VecDeriv InVecDeriv; - typedef linearalgebra::EigenSparseMatrix SparseMatrixEigen; - typedef linearalgebra::EigenSparseMatrix SparseKMatrixEigen; - enum {Nin = In::deriv_total_size, Nout = Out::deriv_total_size }; + SOFA_CLASS(SOFA_TEMPLATE2(DistanceFromTargetMapping,TIn,TOut), SOFA_TEMPLATE3(AssembledNonLinearMapping,TIn,TOut, true)); + + using In = TIn; + using Out = TOut; + using Real = Real_t; + static constexpr auto Nin = In::deriv_total_size; + static constexpr auto Nout = Out::deriv_total_size; + + using InCoord = Coord_t; + using InVecCoord = VecCoord_t; + typedef type::Vec Direction; SOFA_ATTRIBUTE_DEPRECATED__RENAME_DATA_IN_MAPPING_NONLINEAR() @@ -112,21 +103,8 @@ class DistanceFromTargetMapping : public core::Mapping, public BaseDi void init() override; - void apply(const core::MechanicalParams *mparams, Data& out, const Data& in) override; - - void applyJ(const core::MechanicalParams *mparams, Data& out, const Data& in) override; - - void applyJT(const core::MechanicalParams *mparams, Data& out, const Data& in) override; - - void applyJT(const core::ConstraintParams *cparams, Data& out, const Data& in) override; - - void applyDJT(const core::MechanicalParams* mparams, core::MultiVecDerivId parentForce, core::ConstMultiVecDerivId childForce ) override; - - const sofa::linearalgebra::BaseMatrix* getJ() override; - virtual const type::vector* getJs() override; - + void apply(const core::MechanicalParams *mparams, DataVecCoord_t& out, const DataVecCoord_t& in) override; void updateK( const core::MechanicalParams* mparams, core::ConstMultiVecDerivId childForce ) override; - const linearalgebra::BaseMatrix* getK() override; void buildGeometricStiffnessMatrix(sofa::core::GeometricStiffnessMatrix* matrices) override; void draw(const core::visual::VisualParams* vparams) override; @@ -137,9 +115,11 @@ class DistanceFromTargetMapping : public core::Mapping, public BaseDi DistanceFromTargetMapping(); ~DistanceFromTargetMapping() override; - SparseMatrixEigen jacobian; ///< Jacobian of the mapping - type::vector baseMatrices; ///< Jacobian of the mapping, in a vector - SparseKMatrixEigen K; ///< Assembled geometric stiffness matrix + void matrixFreeApplyDJT(const core::MechanicalParams* mparams, Real kFactor, + Data >& parentForce, + const Data >& parentDisplacement, + const Data >& childForce) override; + type::vector directions; ///< Unit vectors in the directions of the lines type::vector< Real > invlengths; ///< inverse of current distances. Null represents the infinity (null distance) diff --git a/Sofa/Component/Mapping/NonLinear/src/sofa/component/mapping/nonlinear/DistanceFromTargetMapping.inl b/Sofa/Component/Mapping/NonLinear/src/sofa/component/mapping/nonlinear/DistanceFromTargetMapping.inl index 79225a2b2c9..1480c7b8301 100644 --- a/Sofa/Component/Mapping/NonLinear/src/sofa/component/mapping/nonlinear/DistanceFromTargetMapping.inl +++ b/Sofa/Component/Mapping/NonLinear/src/sofa/component/mapping/nonlinear/DistanceFromTargetMapping.inl @@ -22,6 +22,7 @@ #pragma once #include +#include #include #include #include @@ -35,12 +36,11 @@ namespace sofa::component::mapping::nonlinear template DistanceFromTargetMapping::DistanceFromTargetMapping() - : Inherit() - , d_indices(initData(&d_indices, "indices", "Indices of the parent points")) + : d_indices(initData(&d_indices, "indices", "Indices of the parent points")) , d_targetPositions(initData(&d_targetPositions, "targetPositions", "Positions to compute the distances from")) , d_restDistances(initData(&d_restDistances, "restLengths", "Rest lengths of the connections")) , d_showObjectScale(initData(&d_showObjectScale, 0.f, "showObjectScale", "Scale for object display")) - , d_color(initData(&d_color, sofa::type::RGBAColor(1,1,0,1), "showColor", "Color for object display. (default=[1.0,1.0,0.0,1.0])")) + , d_color(initData(&d_color, sofa::type::RGBAColor::yellow(), "showColor", "Color for object display.")) { f_indices.setOriginalData(&d_indices); f_targetPositions.setOriginalData(&d_targetPositions); @@ -122,7 +122,7 @@ void DistanceFromTargetMapping::init() baseMatrices.resize( 1 ); baseMatrices[0] = &jacobian; - this->Inherit::init(); // applies the mapping, so after the Data init + Inherit1::init(); // applies the mapping, so after the Data init } template @@ -132,15 +132,16 @@ void DistanceFromTargetMapping::computeCoordPositionDifference( Direc } template -void DistanceFromTargetMapping::apply(const core::MechanicalParams * /*mparams*/ , Data& dOut, const Data& dIn) +void DistanceFromTargetMapping::apply(const core::MechanicalParams * /*mparams*/ , DataVecCoord_t& dOut, const DataVecCoord_t& dIn) { - helper::WriteAccessor< Data > out = dOut; - helper::ReadAccessor< Data > in = dIn; - helper::WriteAccessor > > restDistances(d_restDistances); - const helper::ReadAccessor< Data > > indices(d_indices); - helper::ReadAccessor< Data > targetPositions(d_targetPositions); + helper::WriteAccessor out(dOut); + helper::WriteAccessor restDistances(d_restDistances); + + const helper::ReadAccessor in(dIn); + const helper::ReadAccessor indices(d_indices); + const helper::ReadAccessor targetPositions(d_targetPositions); - jacobian.resizeBlocks(out.size(),in.size()); + this->jacobian.resizeBlocks(out.size(),in.size()); directions.resize(out.size()); invlengths.resize(out.size()); @@ -180,59 +181,26 @@ void DistanceFromTargetMapping::apply(const core::MechanicalParams * jacobian.compress(); } - - -template -void DistanceFromTargetMapping::applyJ(const core::MechanicalParams * /*mparams*/ , Data& dOut, const Data& dIn) -{ - if( jacobian.rowSize() > 0 ) - { - auto dOutWa = sofa::helper::getWriteOnlyAccessor(dOut); - auto dInRa = sofa::helper::getReadAccessor(dIn); - jacobian.mult(dOutWa.wref(),dInRa.ref()); - } - -} - -template -void DistanceFromTargetMapping::applyJT(const core::MechanicalParams * /*mparams*/ , Data& dIn, const Data& dOut) -{ - if( jacobian.rowSize() > 0 ) - { - auto dOutRa = sofa::helper::getReadAccessor(dOut); - auto dInWa = sofa::helper::getWriteOnlyAccessor(dIn); - jacobian.addMultTranspose(dInWa.wref(),dOutRa.ref()); - } -} - -template -void DistanceFromTargetMapping::applyJT(const core::ConstraintParams* cparams, Data& out, const Data& in) -{ - SOFA_UNUSED(cparams); - auto childMatRa = sofa::helper::getReadAccessor(in); - auto parentMatWa = sofa::helper::getWriteAccessor(out); - addMultTransposeEigen(parentMatWa.wref(), jacobian.compressedMatrix, childMatRa.ref()); -} - - template -void DistanceFromTargetMapping::applyDJT(const core::MechanicalParams* mparams, core::MultiVecDerivId parentDfId, core::ConstMultiVecDerivId ) +void DistanceFromTargetMapping::matrixFreeApplyDJT( + const core::MechanicalParams* mparams, Real kFactor, + Data>& parentForce, + const Data>& parentDisplacement, + const Data>& childForce) { - const unsigned geometricStiffness = d_geometricStiffness.getValue().getSelectedId(); - if( !geometricStiffness ) return; + const unsigned geometricStiffness = this->d_geometricStiffness.getValue().getSelectedId(); - helper::WriteAccessor > parentForce (*parentDfId[this->fromModel.get()].write()); - helper::ReadAccessor > parentDisplacement (*mparams->readDx(this->fromModel.get())); // parent displacement - const SReal kfactor = mparams->kFactor(); - helper::ReadAccessor > childForce (*mparams->readF(this->toModel.get())); - const helper::ReadAccessor< Data > > indices(d_indices); + helper::WriteAccessor parentForceAccessor(parentForce); + const helper::ReadAccessor parentDisplacementAccessor(parentDisplacement); + const helper::ReadAccessor childForceAccessor (childForce); + const helper::ReadAccessor indices(d_indices); - for(unsigned i=0; i0) can lead to negative eigen values in geometric stiffness // this results in an undefinite implicit matrix that causes instabilities // if stabilized GS (geometricStiffness==2) -> keep only force in extension - if( childForce[i][0] < 0 || geometricStiffness==1 ) + if( childForceAccessor[i][0] < 0 || geometricStiffness==1 ) { sofa::type::Mat b; // = (I - uu^T) for(unsigned j=0; j::applyDJT(const core::MechanicalParams } } // (I - uu^T)*f/l*kfactor -- do not forget kfactor ! - b *= (Real)(childForce[i][0] * invlengths[i] * kfactor); + b *= (Real)(childForceAccessor[i][0] * invlengths[i] * kFactor); // note that computing a block is not efficient here, but it would // makes sense for storing a stiffness matrix - InDeriv dx = parentDisplacement[indices[i]]; - InDeriv df; + const auto& dx = parentDisplacementAccessor[indices[i]]; + Deriv_t df; for(unsigned j=0; j::applyDJT(const core::MechanicalParams df[j]+=b[j][k]*dx[k]; } } - // InDeriv df = b*dx; - parentForce[indices[i]] += df; + // Deriv_t df = b*dx; + parentForceAccessor[indices[i]] += df; } } } - - - -template -const sofa::linearalgebra::BaseMatrix* DistanceFromTargetMapping::getJ() -{ - return &jacobian; -} - -template -const type::vector* DistanceFromTargetMapping::getJs() -{ - return &baseMatrices; -} - -template -const linearalgebra::BaseMatrix* DistanceFromTargetMapping::getK() -{ - return &K; -} - template void DistanceFromTargetMapping::buildGeometricStiffnessMatrix( sofa::core::GeometricStiffnessMatrix* matrices) { - const unsigned geometricStiffness = d_geometricStiffness.getValue().getSelectedId(); + const unsigned geometricStiffness = this->d_geometricStiffness.getValue().getSelectedId(); if( !geometricStiffness ) { return; @@ -299,7 +246,7 @@ void DistanceFromTargetMapping::buildGeometricStiffnessMatrix( for(sofa::Size i=0; i0) can lead to negative eigen values in geometric stiffness // this results in an undefinite implicit matrix that causes instabilities @@ -327,12 +274,12 @@ template void DistanceFromTargetMapping::updateK( const core::MechanicalParams* mparams, core::ConstMultiVecDerivId childForceId ) { SOFA_UNUSED(mparams); - const unsigned geometricStiffness = d_geometricStiffness.getValue().getSelectedId(); + const unsigned geometricStiffness = this->d_geometricStiffness.getValue().getSelectedId(); if( !geometricStiffness ) { K.resize(0,0); return; } - helper::ReadAccessor > childForce( *childForceId[this->toModel.get()].read() ); - const helper::ReadAccessor< Data > > indices(d_indices); - helper::ReadAccessor > in (*this->fromModel->read(core::ConstVecCoordId::position())); + helper::ReadAccessor childForce( *childForceId[this->toModel.get()].read() ); + const helper::ReadAccessor indices(d_indices); + helper::ReadAccessor in (*this->fromModel->read(core::ConstVecCoordId::position())); K.resizeBlocks(in.size(),in.size()); for(size_t i=0; i Date: Wed, 18 Sep 2024 15:24:04 +0200 Subject: [PATCH 03/16] refactor DistanceMapping --- .../nonlinear/AssembledNonLinearMapping.h | 2 +- .../mapping/nonlinear/DistanceMapping.h | 66 ++----- .../mapping/nonlinear/DistanceMapping.inl | 174 ++++++------------ 3 files changed, 78 insertions(+), 164 deletions(-) diff --git a/Sofa/Component/Mapping/NonLinear/src/sofa/component/mapping/nonlinear/AssembledNonLinearMapping.h b/Sofa/Component/Mapping/NonLinear/src/sofa/component/mapping/nonlinear/AssembledNonLinearMapping.h index 3c7845fa129..86990facad1 100644 --- a/Sofa/Component/Mapping/NonLinear/src/sofa/component/mapping/nonlinear/AssembledNonLinearMapping.h +++ b/Sofa/Component/Mapping/NonLinear/src/sofa/component/mapping/nonlinear/AssembledNonLinearMapping.h @@ -79,7 +79,7 @@ class AssembledNonLinearMapping : public core::Mapping, public NonLin struct JacobianEntry { sofa::Index vertexId; - typename In::Coord jacobianValue; + typename In::CPos jacobianValue; bool operator<(const JacobianEntry& other) const { return vertexId < other.vertexId;} }; }; diff --git a/Sofa/Component/Mapping/NonLinear/src/sofa/component/mapping/nonlinear/DistanceMapping.h b/Sofa/Component/Mapping/NonLinear/src/sofa/component/mapping/nonlinear/DistanceMapping.h index 07ea6194f49..caee530a3cb 100644 --- a/Sofa/Component/Mapping/NonLinear/src/sofa/component/mapping/nonlinear/DistanceMapping.h +++ b/Sofa/Component/Mapping/NonLinear/src/sofa/component/mapping/nonlinear/DistanceMapping.h @@ -23,7 +23,7 @@ #include -#include +#include #include #include #include @@ -46,26 +46,18 @@ namespace sofa::component::mapping::nonlinear * @author Francois Faure */ template -class DistanceMapping : public core::Mapping, public StabilizedNonLinearMappingData +class DistanceMapping : public AssembledNonLinearMapping { public: - SOFA_CLASS(SOFA_TEMPLATE2(DistanceMapping,TIn,TOut), SOFA_TEMPLATE2(core::Mapping,TIn,TOut)); + SOFA_CLASS(SOFA_TEMPLATE2(DistanceMapping,TIn,TOut), SOFA_TEMPLATE3(AssembledNonLinearMapping,TIn,TOut,true)); typedef TIn In; typedef TOut Out; - typedef typename Out::VecCoord OutVecCoord; - typedef typename Out::VecDeriv OutVecDeriv; - typedef typename Out::Deriv OutDeriv; - typedef typename Out::MatrixDeriv OutMatrixDeriv; - typedef typename Out::Real Real; - typedef typename In::Deriv InDeriv; - typedef typename In::MatrixDeriv InMatrixDeriv; - typedef typename In::Coord InCoord; - typedef typename In::VecCoord InVecCoord; - typedef typename In::VecDeriv InVecDeriv; - typedef linearalgebra::EigenSparseMatrix SparseMatrixEigen; - typedef linearalgebra::EigenSparseMatrix SparseKMatrixEigen; + + using Real = Real_t; + static constexpr auto Nin = In::deriv_total_size; + typedef sofa::core::topology::BaseMeshTopology::SeqEdges SeqEdges; typedef type::Vec Direction; @@ -84,25 +76,10 @@ class DistanceMapping : public core::Mapping, public StabilizedNonLin /// Link to be set to the topology container in the component graph. SingleLink, sofa::core::topology::BaseMeshTopology, BaseLink::FLAG_STOREPATH | BaseLink::FLAG_STRONGLINK> l_topology; - - void init() override; - void apply(const core::MechanicalParams *mparams, Data& out, const Data& in) override; - - void applyJ(const core::MechanicalParams *mparams, Data& out, const Data& in) override; - - void applyJT(const core::MechanicalParams *mparams, Data& out, const Data& in) override; - - void applyJT(const core::ConstraintParams *cparams, Data& out, const Data& in) override; - - void applyDJT(const core::MechanicalParams* mparams, core::MultiVecDerivId parentForce, core::ConstMultiVecDerivId childForce ) override; - - const sofa::linearalgebra::BaseMatrix* getJ() override; - virtual const type::vector* getJs() override; - + void apply(const core::MechanicalParams *mparams, DataVecCoord_t& out, const DataVecCoord_t& in) override; void updateK( const core::MechanicalParams* mparams, core::ConstMultiVecDerivId childForce ) override; - const linearalgebra::BaseMatrix* getK() override; void buildGeometricStiffnessMatrix(sofa::core::GeometricStiffnessMatrix* matrices) override; void computeBBox(const core::ExecParams* params, bool onlyVisible) override; @@ -111,29 +88,18 @@ class DistanceMapping : public core::Mapping, public StabilizedNonLin protected: DistanceMapping(); - SparseMatrixEigen jacobian; ///< Jacobian of the mapping - type::vector baseMatrices; ///< Jacobian of the mapping, in a vector - SparseKMatrixEigen K; ///< Assembled geometric stiffness matrix + void matrixFreeApplyDJT(const core::MechanicalParams* mparams, Real kFactor, + Data >& parentForce, + const Data >& parentDisplacement, + const Data >& childForce) override; + type::vector directions; ///< Unit vectors in the directions of the lines type::vector< Real > invlengths; ///< inverse of current distances. Null represents the infinity (null distance) /// r=b-a only for position (eventual rotation, affine transform... remains null) - void computeCoordPositionDifference( Direction& r, const InCoord& a, const InCoord& b ); - - /** - * @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; - Direction jacobianValue; - bool operator<(const JacobianEntry& other) const { return vertexId < other.vertexId;} - }; + void computeCoordPositionDifference( Direction& r, const Coord_t& a, const Coord_t& b ); + + using Inherit1::JacobianEntry; }; diff --git a/Sofa/Component/Mapping/NonLinear/src/sofa/component/mapping/nonlinear/DistanceMapping.inl b/Sofa/Component/Mapping/NonLinear/src/sofa/component/mapping/nonlinear/DistanceMapping.inl index c4b91bb0a4a..10c682f833f 100644 --- a/Sofa/Component/Mapping/NonLinear/src/sofa/component/mapping/nonlinear/DistanceMapping.inl +++ b/Sofa/Component/Mapping/NonLinear/src/sofa/component/mapping/nonlinear/DistanceMapping.inl @@ -22,6 +22,7 @@ #pragma once #include +#include #include #include #include @@ -69,7 +70,7 @@ void DistanceMapping::init() typename core::behavior::MechanicalState::ReadVecCoord pos = this->getFromModel()->readPositions(); this->getToModel()->resize( links.size() ); - jacobian.resizeBlocks(links.size(),pos.size()); + this->jacobian.resizeBlocks(links.size(),pos.size()); directions.resize(links.size()); invlengths.resize(links.size()); @@ -90,27 +91,27 @@ void DistanceMapping::init() } } - baseMatrices.resize( 1 ); - baseMatrices[0] = &jacobian; + this->baseMatrices.resize( 1 ); + this->baseMatrices[0] = &this->jacobian; this->Inherit1::init(); // applies the mapping, so after the Data init } template -void DistanceMapping::computeCoordPositionDifference( Direction& r, const InCoord& a, const InCoord& b ) +void DistanceMapping::computeCoordPositionDifference( Direction& r, const Coord_t& a, const Coord_t& b ) { r = TIn::getCPos(b) - TIn::getCPos(a); } template -void DistanceMapping::apply(const core::MechanicalParams * /*mparams*/ , Data& dOut, const Data& dIn) +void DistanceMapping::apply(const core::MechanicalParams * /*mparams*/ , DataVecCoord_t& dOut, const DataVecCoord_t& dIn) { - helper::WriteOnlyAccessor< Data > out = dOut; - helper::ReadAccessor< Data > in = dIn; + helper::WriteOnlyAccessor> out(dOut); + helper::ReadAccessor in(dIn); helper::ReadAccessor > > restLengths(d_restLengths); const SeqEdges& links = l_topology->getEdges(); - jacobian.clear(); + this->jacobian.clear(); for (unsigned int i = 0; i < links.size(); ++i) { @@ -138,139 +139,92 @@ void DistanceMapping::apply(const core::MechanicalParams * /*mparams* direction.fill(p); } - sofa::type::fixed_array jacobianEntries {JacobianEntry{link[0], -direction}, JacobianEntry{link[1], direction}}; + sofa::type::fixed_array jacobianEntries { + JacobianEntry{link[0], -direction}, + JacobianEntry{link[1], direction} + }; //invert to insert in increasing column order std::sort(jacobianEntries.begin(), jacobianEntries.end()); - jacobian.beginRow(i); + this->jacobian.beginRow(i); for (const auto& [vertexId, jacobianValue] : jacobianEntries) { for (unsigned k = 0; k < In::spatial_dimensions; ++k) { - jacobian.insertBack(i, vertexId * Nin + k, jacobianValue[k]); + this->jacobian.insertBack(i, vertexId * Nin + k, jacobianValue[k]); } } } - jacobian.compress(); + this->jacobian.compress(); } - template -void DistanceMapping::applyJ(const core::MechanicalParams * /*mparams*/ , Data& out, const Data& in) +void DistanceMapping::matrixFreeApplyDJT( + const core::MechanicalParams* mparams, Real kFactor, + Data>& parentForce, + const Data>& parentDisplacement, + const Data>& childForce) { - if( jacobian.rowSize() ) - { - auto dOutWa = sofa::helper::getWriteOnlyAccessor(out); - auto dInRa = sofa::helper::getReadAccessor(in); - jacobian.mult(dOutWa.wref(),dInRa.ref()); - } -} + SOFA_UNUSED(mparams); + const unsigned geometricStiffness = this->d_geometricStiffness.getValue().getSelectedId(); -template -void DistanceMapping::applyJT(const core::MechanicalParams * /*mparams*/ , Data& out, const Data& in) -{ - if( jacobian.rowSize() ) - { - auto dOutRa = sofa::helper::getReadAccessor(in); - auto dInWa = sofa::helper::getWriteOnlyAccessor(out); - jacobian.addMultTranspose(dInWa.wref(),dOutRa.ref()); - } -} + helper::WriteAccessor parentForceAccessor(parentForce); + helper::ReadAccessor parentDisplacementAccessor(parentDisplacement); + helper::ReadAccessor childForceAccessor(childForce); -template -void DistanceMapping::applyDJT(const core::MechanicalParams* mparams, core::MultiVecDerivId parentDfId, core::ConstMultiVecDerivId) -{ - const unsigned geometricStiffness = d_geometricStiffness.getValue().getSelectedId(); - if( !geometricStiffness ) return; - - helper::WriteAccessor > parentForce (*parentDfId[this->fromModel.get()].write()); - helper::ReadAccessor > parentDisplacement (*mparams->readDx(this->fromModel.get())); // parent displacement - const SReal& kfactor = mparams->kFactor(); - helper::ReadAccessor > childForce (*mparams->readF(this->toModel.get())); + const SeqEdges& links = l_topology->getEdges(); - if( K.compressedMatrix.nonZeros() ) - { - K.addMult( parentForce.wref(), parentDisplacement.ref(), (typename In::Real)kfactor ); - } - else + for (unsigned i = 0; i < links.size(); ++i) { - const SeqEdges& links = l_topology->getEdges(); - - for(unsigned i=0; i0) can lead to negative eigen values in geometric stiffness + // this results in an undefinite implicit matrix that causes instabilities + // if stabilized GS (geometricStiffness==2) -> keep only force in extension + if( childForceAccessor[i][0] < 0 || geometricStiffness==1 ) { - // force in compression (>0) can lead to negative eigen values in geometric stiffness - // this results in an undefinite implicit matrix that causes instabilities - // if stabilized GS (geometricStiffness==2) -> keep only force in extension - if( childForce[i][0] < 0 || geometricStiffness==1 ) + sofa::type::Mat b; // = (I - uu^T) + for(unsigned j=0; j b; // = (I - uu^T) - for(unsigned j=0; j(1) * ( j==k ) - directions[i][j]*directions[i][k]; - } + b[j][k] = static_cast(1) * ( j==k ) - directions[i][j]*directions[i][k]; } - // (I - uu^T)*f/l*kfactor -- do not forget kfactor ! - b *= (Real)(childForce[i][0] * invlengths[i] * kfactor); - // note that computing a block is not efficient here, but it - // would make sense for storing a stiffness matrix - - InDeriv dx = parentDisplacement[links[i][1]] - parentDisplacement[links[i][0]]; - InDeriv df; - for(unsigned j=0; j dx = parentDisplacementAccessor[links[i][1]] - parentDisplacementAccessor[links[i][0]]; + Deriv_t df; + for(unsigned j=0; j -void DistanceMapping::applyJT(const core::ConstraintParams* cparams, Data& out, const Data& in) -{ - SOFA_UNUSED(cparams); - auto childMatRa = sofa::helper::getReadAccessor(in); - auto parentMatWa = sofa::helper::getWriteAccessor(out); - addMultTransposeEigen(parentMatWa.wref(), jacobian.compressedMatrix, childMatRa.ref()); -} - - -template -const sofa::linearalgebra::BaseMatrix* DistanceMapping::getJ() -{ - return &jacobian; -} - -template -const type::vector* DistanceMapping::getJs() -{ - return &baseMatrices; -} - - template void DistanceMapping::updateK(const core::MechanicalParams *mparams, core::ConstMultiVecDerivId childForceId ) { 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(); + if( !geometricStiffness ) { this->K.resize(0,0); return; } - helper::ReadAccessor > childForce( *childForceId[this->toModel.get()].read() ); + const helper::ReadAccessor childForce( *childForceId[this->toModel.get()].read() ); const SeqEdges& links = l_topology->getEdges(); unsigned int size = this->fromModel->getSize(); - K.resizeBlocks(size,size); + this->K.resizeBlocks(size,size); for(size_t i=0; i0) can lead to negative eigen values in geometric stiffness @@ -290,26 +244,20 @@ void DistanceMapping::updateK(const core::MechanicalParams *mparams, b *= childForce[i][0] * invlengths[i]; // (I - uu^T)*f/l // Note that 'links' is not sorted so the matrix can not be filled-up in order - K.addBlock(links[i][0],links[i][0],b); - K.addBlock(links[i][0],links[i][1],-b); - K.addBlock(links[i][1],links[i][0],-b); - K.addBlock(links[i][1],links[i][1],b); + this->K.addBlock(links[i][0],links[i][0],b); + this->K.addBlock(links[i][0],links[i][1],-b); + this->K.addBlock(links[i][1],links[i][0],-b); + this->K.addBlock(links[i][1],links[i][1],b); } } - K.compress(); -} - -template -const linearalgebra::BaseMatrix* DistanceMapping::getK() -{ - return &K; + this->K.compress(); } template void DistanceMapping::buildGeometricStiffnessMatrix( sofa::core::GeometricStiffnessMatrix* matrices) { - const unsigned& geometricStiffness = d_geometricStiffness.getValue().getSelectedId(); + const unsigned& geometricStiffness = this->d_geometricStiffness.getValue().getSelectedId(); if( !geometricStiffness ) { return; @@ -321,7 +269,7 @@ void DistanceMapping::buildGeometricStiffnessMatrix( for(sofa::Size i=0; i force_i = childForce[i]; // force in compression (>0) can lead to negative eigen values in geometric stiffness // this results in an undefinite implicit matrix that causes instabilities From 0cbd6a18343cee0967dbe7c710edcb0e53bf46bc Mon Sep 17 00:00:00 2001 From: Alex Bilger Date: Wed, 18 Sep 2024 15:33:53 +0200 Subject: [PATCH 04/16] refactor SquareDistanceMapping --- .../mapping/nonlinear/SquareDistanceMapping.h | 68 ++----- .../nonlinear/SquareDistanceMapping.inl | 169 ++++++------------ 2 files changed, 74 insertions(+), 163 deletions(-) diff --git a/Sofa/Component/Mapping/NonLinear/src/sofa/component/mapping/nonlinear/SquareDistanceMapping.h b/Sofa/Component/Mapping/NonLinear/src/sofa/component/mapping/nonlinear/SquareDistanceMapping.h index afca0b322b5..491cdbcb699 100644 --- a/Sofa/Component/Mapping/NonLinear/src/sofa/component/mapping/nonlinear/SquareDistanceMapping.h +++ b/Sofa/Component/Mapping/NonLinear/src/sofa/component/mapping/nonlinear/SquareDistanceMapping.h @@ -23,8 +23,7 @@ #include -#include -#include +#include #include #include #include @@ -53,38 +52,21 @@ namespace sofa::component::mapping::nonlinear // If the rest lengths are not defined, they are set using the initial values. // If computeDistance is set to true, the rest lengths are set to 0. template -class SquareDistanceMapping : public core::Mapping, public NonLinearMappingData +class SquareDistanceMapping : public AssembledNonLinearMapping { public: - SOFA_CLASS(SOFA_TEMPLATE2(SquareDistanceMapping,TIn,TOut), SOFA_TEMPLATE2(core::Mapping,TIn,TOut)); - - typedef core::Mapping Inherit; - typedef TIn In; - typedef TOut Out; - typedef typename Out::VecCoord OutVecCoord; - typedef typename Out::VecDeriv OutVecDeriv; - typedef typename Out::Coord OutCoord; - typedef typename Out::Deriv OutDeriv; - typedef typename Out::MatrixDeriv OutMatrixDeriv; - typedef typename Out::Real Real; - typedef typename In::Deriv InDeriv; - typedef typename In::MatrixDeriv InMatrixDeriv; - typedef typename In::Coord InCoord; - typedef typename In::VecCoord InVecCoord; - typedef typename In::VecDeriv InVecDeriv; - typedef linearalgebra::EigenSparseMatrix SparseMatrixEigen; - typedef linearalgebra::EigenSparseMatrix SparseKMatrixEigen; - typedef Data InDataVecCoord; - typedef Data InDataVecDeriv; - typedef Data InDataMatrixDeriv; - typedef Data OutDataVecCoord; - typedef Data OutDataVecDeriv; - typedef Data OutDataMatrixDeriv; - enum {Nin = In::deriv_total_size, Nout = Out::deriv_total_size }; + SOFA_CLASS(SOFA_TEMPLATE2(SquareDistanceMapping,TIn,TOut), SOFA_TEMPLATE3(AssembledNonLinearMapping,TIn,TOut,true)); + + using In = TIn; + using Out = TOut; + + using Real = Real_t; + + static constexpr auto Nin = In::deriv_total_size; + typedef sofa::core::topology::BaseMeshTopology::SeqEdges SeqEdges; typedef type::Vec Direction; - Data d_showObjectScale; ///< Scale for object display Data d_color; ///< Color for object display. (default=[1.0,1.0,0.0,1.0]) @@ -93,37 +75,23 @@ class SquareDistanceMapping : public core::Mapping, public NonLinearM void init() override; - using Inherit::apply; - - void apply(const core::MechanicalParams *mparams, Data& out, const Data& in) override; - - void applyJ(const core::MechanicalParams *mparams, Data& out, const Data& in) override; - - void applyJT(const core::MechanicalParams *mparams, Data& out, const Data& in) override; - - void applyJT(const core::ConstraintParams *cparams, Data& out, const Data& in) override; - - void applyDJT(const core::MechanicalParams* mparams, core::MultiVecDerivId parentForce, core::ConstMultiVecDerivId childForce ) override; - - const sofa::linearalgebra::BaseMatrix* getJ() override; - virtual const type::vector* getJs() override; - + void apply(const core::MechanicalParams *mparams, DataVecCoord_t& out, const DataVecCoord_t& in) override; void updateK( const core::MechanicalParams* mparams, core::ConstMultiVecDerivId childForce ) override; - const linearalgebra::BaseMatrix* getK() override; void buildGeometricStiffnessMatrix(sofa::core::GeometricStiffnessMatrix* matrices) override; void draw(const core::visual::VisualParams* vparams) override; protected: SquareDistanceMapping(); - virtual ~SquareDistanceMapping(); + ~SquareDistanceMapping() override; - SparseMatrixEigen jacobian; ///< Jacobian of the mapping - type::vector baseMatrices; ///< Jacobian of the mapping, in a vector - SparseKMatrixEigen K; ///< Assembled geometric stiffness matrix + void matrixFreeApplyDJT(const core::MechanicalParams* mparams, Real kFactor, + Data >& parentForce, + const Data >& parentDisplacement, + const Data >& childForce) override; /// r=b-a only for position (eventual rotation, affine transform... remains null) - void computeCoordPositionDifference( Direction& r, const InCoord& a, const InCoord& b ); + void computeCoordPositionDifference( Direction& r, const Coord_t& a, const Coord_t& b ); }; diff --git a/Sofa/Component/Mapping/NonLinear/src/sofa/component/mapping/nonlinear/SquareDistanceMapping.inl b/Sofa/Component/Mapping/NonLinear/src/sofa/component/mapping/nonlinear/SquareDistanceMapping.inl index ae9a7e81e51..a88a9e02647 100644 --- a/Sofa/Component/Mapping/NonLinear/src/sofa/component/mapping/nonlinear/SquareDistanceMapping.inl +++ b/Sofa/Component/Mapping/NonLinear/src/sofa/component/mapping/nonlinear/SquareDistanceMapping.inl @@ -22,6 +22,7 @@ #pragma once #include +#include #include #include #include @@ -36,8 +37,7 @@ namespace sofa::component::mapping::nonlinear template SquareDistanceMapping::SquareDistanceMapping() - : Inherit() - , d_showObjectScale(initData(&d_showObjectScale, Real(0), "showObjectScale", "Scale for object display")) + : d_showObjectScale(initData(&d_showObjectScale, Real(0), "showObjectScale", "Scale for object display")) , d_color(initData(&d_color, sofa::type::RGBAColor(1,1,0,1), "showColor", "Color for object display. (default=[1.0,1.0,0.0,1.0])")) , l_topology(initLink("topology", "link to the topology container")) { @@ -45,9 +45,7 @@ SquareDistanceMapping::SquareDistanceMapping() template SquareDistanceMapping::~SquareDistanceMapping() -{ -} - += default; template void SquareDistanceMapping::init() @@ -71,34 +69,37 @@ void SquareDistanceMapping::init() this->getToModel()->resize( links.size() ); - baseMatrices.resize( 1 ); - baseMatrices[0] = &jacobian; + this->baseMatrices.resize( 1 ); + this->baseMatrices[0] = &this->jacobian; - this->Inherit::init(); // applies the mapping, so after the Data init + this->Inherit1::init(); } template -void SquareDistanceMapping::computeCoordPositionDifference( Direction& r, const InCoord& a, const InCoord& b ) +void SquareDistanceMapping::computeCoordPositionDifference( Direction& r, const Coord_t& a, const Coord_t& b ) { r = TIn::getCPos(b)-TIn::getCPos(a); } template -void SquareDistanceMapping::apply(const core::MechanicalParams * /*mparams*/ , Data& dOut, const Data& dIn) +void SquareDistanceMapping::apply( + const core::MechanicalParams * /*mparams*/, + DataVecCoord_t& dOut, + const DataVecCoord_t& dIn) { - helper::WriteOnlyAccessor< Data > out = dOut; - helper::ReadAccessor< Data > in = dIn; + helper::WriteOnlyAccessor< DataVecCoord_t > out = dOut; + const helper::ReadAccessor in (dIn); const SeqEdges& links = l_topology->getEdges(); - jacobian.resizeBlocks(out.size(),in.size()); + this->jacobian.resizeBlocks(out.size(),in.size()); Direction gap; for(unsigned i=0; i& p0 = in[links[i][0]]; + const Coord_t& p1 = in[links[i][1]]; // gap = in[links[i][1]] - in[links[i][0]] (only for position) computeCoordPositionDifference( gap, p0, p1 ); @@ -120,126 +121,74 @@ void SquareDistanceMapping::apply(const core::MechanicalParams * /*mp // } - jacobian.beginRow(i); + this->jacobian.beginRow(i); if( links[i][1]jacobian.insertBack( i, links[i][1]*Nin+k, gap[k] ); for(unsigned k=0; kjacobian.insertBack( i, links[i][0]*Nin+k, -gap[k] ); } else { for(unsigned k=0; kjacobian.insertBack( i, links[i][0]*Nin+k, -gap[k] ); for(unsigned k=0; kjacobian.insertBack( i, links[i][1]*Nin+k, gap[k] ); } } - jacobian.compress(); -} - - -template -void SquareDistanceMapping::applyJ(const core::MechanicalParams * /*mparams*/ , Data& dOut, const Data& dIn) -{ - if( jacobian.rowSize() ) - { - auto dOutWa = sofa::helper::getWriteOnlyAccessor(dOut); - auto dInRa = sofa::helper::getReadAccessor(dIn); - jacobian.mult(dOutWa.wref(),dInRa.ref()); - } + this->jacobian.compress(); } template -void SquareDistanceMapping::applyJT(const core::MechanicalParams * /*mparams*/ , Data& dIn, const Data& dOut) +void SquareDistanceMapping::matrixFreeApplyDJT( + const core::MechanicalParams* mparams, Real kFactor, + Data>& parentForce, + const Data>& parentDisplacement, + const Data>& childForce) { - if( jacobian.rowSize() ) - { - auto dOutRa = sofa::helper::getReadAccessor(dOut); - auto dInWa = sofa::helper::getWriteOnlyAccessor(dIn); - jacobian.addMultTranspose(dInWa.wref(),dOutRa.ref()); - } -} + SOFA_UNUSED(mparams); + const unsigned& geometricStiffness = this->d_geometricStiffness.getValue().getSelectedId(); -template -void SquareDistanceMapping::applyDJT(const core::MechanicalParams* mparams, core::MultiVecDerivId parentDfId, core::ConstMultiVecDerivId ) -{ - const unsigned& geometricStiffness = d_geometricStiffness.getValue().getSelectedId(); - if( !geometricStiffness ) return; + helper::WriteAccessor parentForceAccessor(parentForce); + helper::ReadAccessor parentDisplacementAccessor(parentDisplacement); + helper::ReadAccessor childForceAccessor(childForce); - helper::WriteAccessor > parentForce (*parentDfId[this->fromModel.get()].write()); - helper::ReadAccessor > parentDisplacement (*mparams->readDx(this->fromModel.get())); // parent displacement - const SReal& kfactor = mparams->kFactor(); - helper::ReadAccessor > childForce (*mparams->readF(this->toModel.get())); + const SeqEdges& links = l_topology->getEdges(); - if( K.compressedMatrix.nonZeros() ) + for (unsigned i = 0; i < links.size(); i++) { - K.addMult( parentForce.wref(), parentDisplacement.ref(), (typename In::Real)kfactor ); - } - else - { - const SeqEdges& links = l_topology->getEdges(); - - for(unsigned i=0; i0) can lead to negative eigen values in geometric stiffness + // this results in an undefinite implicit matrix that causes instabilities + // if stabilized GS (geometricStiffness==2) -> keep only force in extension + if( childForceAccessor[i][0] < 0 || geometricStiffness==1 ) { - // force in compression (>0) can lead to negative eigen values in geometric stiffness - // this results in an undefinite implicit matrix that causes instabilities - // if stabilized GS (geometricStiffness==2) -> keep only force in extension - if( childForce[i][0] < 0 || geometricStiffness==1 ) - { - - SReal tmp = 2*childForce[i][0]*kfactor; + SReal tmp = 2 * childForceAccessor[i][0] * kFactor; - typename In::DPos df = tmp * ( - In::getDPos(parentDisplacement[links[i][0]]) - - In::getDPos(parentDisplacement[links[i][1]])); - // it is symmetric so -df = (parentDisplacement[links[i][1]]-parentDisplacement[links[i][0]])*tmp; + typename In::DPos df = tmp * ( + In::getDPos(parentDisplacementAccessor[links[i][0]]) - + In::getDPos(parentDisplacementAccessor[links[i][1]])); + // it is symmetric so -df = (parentDisplacement[links[i][1]]-parentDisplacement[links[i][0]])*tmp; - In::setDPos(parentForce[links[i][0]], In::getDPos(parentForce[links[i][0]]) + df); - In::setDPos(parentForce[links[i][1]], In::getDPos(parentForce[links[i][1]]) - df); - } + In::setDPos(parentForceAccessor[links[i][0]], In::getDPos(parentForceAccessor[links[i][0]]) + df); + In::setDPos(parentForceAccessor[links[i][1]], In::getDPos(parentForceAccessor[links[i][1]]) - df); } } } -template -void SquareDistanceMapping::applyJT(const core::ConstraintParams* cparams, Data& out, const Data& in) -{ - SOFA_UNUSED(cparams); - auto childMatRa = sofa::helper::getReadAccessor(in); - auto parentMatWa = sofa::helper::getWriteAccessor(out); - addMultTransposeEigen(parentMatWa.wref(), jacobian.compressedMatrix, childMatRa.ref()); -} - - -template -const sofa::linearalgebra::BaseMatrix* SquareDistanceMapping::getJ() -{ - return &jacobian; -} - -template -const type::vector* SquareDistanceMapping::getJs() -{ - return &baseMatrices; -} - - - template void SquareDistanceMapping::updateK(const core::MechanicalParams *mparams, core::ConstMultiVecDerivId childForceId ) { 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(); + if( !geometricStiffness ) { this->K.resize(0,0); return; } - helper::ReadAccessor > childForce( *childForceId[this->toModel.get()].read() ); + helper::ReadAccessor childForce( *childForceId[this->toModel.get()].read() ); const SeqEdges& links = l_topology->getEdges(); unsigned int size = this->fromModel->getSize(); - K.resizeBlocks(size,size); + this->K.resizeBlocks(size,size); for(size_t i=0; i0) can lead to negative eigen values in geometric stiffness @@ -251,27 +200,21 @@ void SquareDistanceMapping::updateK(const core::MechanicalParams *mpa for(unsigned k=0; kK.add( links[i][0]*Nin+k, links[i][0]*Nin+k, tmp ); + this->K.add( links[i][0]*Nin+k, links[i][1]*Nin+k, -tmp ); + this->K.add( links[i][1]*Nin+k, links[i][1]*Nin+k, tmp ); + this->K.add( links[i][1]*Nin+k, links[i][0]*Nin+k, -tmp ); } } } - K.compress(); -} - -template -const linearalgebra::BaseMatrix* SquareDistanceMapping::getK() -{ - return &K; + this->K.compress(); } template void SquareDistanceMapping::buildGeometricStiffnessMatrix( sofa::core::GeometricStiffnessMatrix* matrices) { - const unsigned geometricStiffness = d_geometricStiffness.getValue().getSelectedId(); + const unsigned geometricStiffness = this->d_geometricStiffness.getValue().getSelectedId(); if( !geometricStiffness ) { return; @@ -283,7 +226,7 @@ void SquareDistanceMapping::buildGeometricStiffnessMatrix( for(sofa::Size i=0; i& force_i = childForce[i]; const sofa::topology::Edge link = links[i]; // force in compression (>0) can lead to negative eigen values in geometric stiffness From 323ad230e8220a5c9714494233f04bf07082c9c7 Mon Sep 17 00:00:00 2001 From: Alex Bilger Date: Wed, 18 Sep 2024 15:55:21 +0200 Subject: [PATCH 05/16] refactor SquareMapping --- .../mapping/nonlinear/SquareMapping.h | 63 ++------ .../mapping/nonlinear/SquareMapping.inl | 144 +++++------------- .../NonLinear/tests/SquareMapping_test.cpp | 10 -- 3 files changed, 48 insertions(+), 169 deletions(-) diff --git a/Sofa/Component/Mapping/NonLinear/src/sofa/component/mapping/nonlinear/SquareMapping.h b/Sofa/Component/Mapping/NonLinear/src/sofa/component/mapping/nonlinear/SquareMapping.h index 161037cbb24..40e89b9a14a 100644 --- a/Sofa/Component/Mapping/NonLinear/src/sofa/component/mapping/nonlinear/SquareMapping.h +++ b/Sofa/Component/Mapping/NonLinear/src/sofa/component/mapping/nonlinear/SquareMapping.h @@ -23,7 +23,7 @@ #include -#include +#include #include #include @@ -38,65 +38,30 @@ namespace sofa::component::mapping::nonlinear */ template -class SquareMapping : public core::Mapping, public NonLinearMappingData +class SquareMapping : public AssembledNonLinearMapping { public: - SOFA_CLASS(SOFA_TEMPLATE2(SquareMapping,TIn,TOut), SOFA_TEMPLATE2(core::Mapping,TIn,TOut)); - - typedef core::Mapping Inherit; - typedef TIn In; - typedef TOut Out; - typedef typename Out::VecCoord OutVecCoord; - typedef typename Out::VecDeriv OutVecDeriv; - typedef typename Out::Coord OutCoord; - typedef typename Out::Deriv OutDeriv; - typedef typename Out::MatrixDeriv OutMatrixDeriv; - typedef typename Out::Real Real; - typedef typename In::Deriv InDeriv; - typedef typename In::MatrixDeriv InMatrixDeriv; - typedef typename In::Coord InCoord; - typedef typename In::VecCoord InVecCoord; - typedef typename In::VecDeriv InVecDeriv; - typedef linearalgebra::EigenSparseMatrix SparseMatrixEigen; - typedef linearalgebra::EigenSparseMatrix SparseKMatrixEigen; - typedef Data InDataVecCoord; - typedef Data InDataVecDeriv; - typedef Data InDataMatrixDeriv; - typedef Data OutDataVecCoord; - typedef Data OutDataVecDeriv; - typedef Data OutDataMatrixDeriv; - typedef type::Vec Direction; + SOFA_CLASS(SOFA_TEMPLATE2(SquareMapping,TIn,TOut), SOFA_TEMPLATE3(AssembledNonLinearMapping,TIn,TOut,false)); - void init() override; - - using Inherit::apply; - - void apply(const core::MechanicalParams *mparams, Data& out, const Data& in) override; - - void applyJ(const core::MechanicalParams *mparams, Data& out, const Data& in) override; - - void applyJT(const core::MechanicalParams *mparams, Data& out, const Data& in) override; + using In = TIn; + using Out = TOut; - void applyJT(const core::ConstraintParams *cparams, Data& out, const Data& in) override; + using Real = Real_t; - void applyDJT(const core::MechanicalParams* mparams, core::MultiVecDerivId parentForce, core::ConstMultiVecDerivId childForce ) override; + static constexpr auto Nin = In::deriv_total_size; - const sofa::linearalgebra::BaseMatrix* getJ() override; - virtual const type::vector* getJs() override; + void init() override; + void apply(const core::MechanicalParams *mparams, DataVecCoord_t& out, const DataVecCoord_t& in) override; void updateK( const core::MechanicalParams* mparams, core::ConstMultiVecDerivId childForce ) override; - const linearalgebra::BaseMatrix* getK() override; void buildGeometricStiffnessMatrix(sofa::core::GeometricStiffnessMatrix* matrices) override; - Data < bool > d_useGeometricStiffnessMatrix; ///< If available (cached), the geometric stiffness matrix is used in order to compute the product with the parent displacement. Otherwise, the product is computed directly using the available vectors (matrix-free method). - protected: - SquareMapping(); - ~SquareMapping() override; - SparseMatrixEigen jacobian; ///< Jacobian of the mapping - type::vector baseMatrices; ///< Jacobian of the mapping, in a vector - SparseKMatrixEigen K; ///< Assembled geometric stiffness matrix + void matrixFreeApplyDJT(const core::MechanicalParams* mparams, Real kFactor, + Data >& parentForce, + const Data >& parentDisplacement, + const Data >& childForce) override; }; @@ -104,8 +69,6 @@ class SquareMapping : public core::Mapping, public NonLinearMappingDa #if !defined(SOFA_COMPONENT_MAPPING_SquareMapping_CPP) extern template class SOFA_COMPONENT_MAPPING_NONLINEAR_API SquareMapping< defaulttype::Vec1Types, defaulttype::Vec1Types >; - - #endif } // namespace sofa::component::mapping::nonlinear diff --git a/Sofa/Component/Mapping/NonLinear/src/sofa/component/mapping/nonlinear/SquareMapping.inl b/Sofa/Component/Mapping/NonLinear/src/sofa/component/mapping/nonlinear/SquareMapping.inl index 3f7eb33a74d..9f69a3893cd 100644 --- a/Sofa/Component/Mapping/NonLinear/src/sofa/component/mapping/nonlinear/SquareMapping.inl +++ b/Sofa/Component/Mapping/NonLinear/src/sofa/component/mapping/nonlinear/SquareMapping.inl @@ -22,6 +22,7 @@ #pragma once #include +#include #include #include #include @@ -32,160 +33,85 @@ namespace sofa::component::mapping::nonlinear { -template -SquareMapping::SquareMapping() - : Inherit() - , d_useGeometricStiffnessMatrix(initData(&d_useGeometricStiffnessMatrix, true, "useGeometricStiffnessMatrix", - "If available (cached), the geometric stiffness matrix is used in order to compute the " - "product with the parent displacement. Otherwise, the product is computed directly using " - "the available vectors (matrix-free method).")) -{ -} - -template -SquareMapping::~SquareMapping() -{ -} - - template void SquareMapping::init() { - baseMatrices.resize( 1 ); - baseMatrices[0] = &jacobian; + this->baseMatrices.resize( 1 ); + this->baseMatrices[0] = &this->jacobian; - this->Inherit::init(); + Inherit1::init(); } - template -void SquareMapping::apply(const core::MechanicalParams * /*mparams*/ , Data& dOut, const Data& dIn) +void SquareMapping::apply(const core::MechanicalParams* mparams, + DataVecCoord_t& dOut, const DataVecCoord_t& dIn) { - helper::WriteOnlyAccessor< Data > out = dOut; - helper::ReadAccessor< Data > in = dIn; + helper::WriteOnlyAccessor< DataVecCoord_t > out = dOut; + const helper::ReadAccessor> in = dIn; size_t size = in.size(); this->getToModel()->resize( size ); - jacobian.resizeBlocks( size, size ); - jacobian.reserve( size ); + this->jacobian.resizeBlocks( size, size ); + this->jacobian.reserve( size ); for( unsigned i=0 ; ijacobian.beginRow(i); + this->jacobian.insertBack( i, i, 2.0*x ); } - jacobian.compress(); -} - - -template -void SquareMapping::applyJ(const core::MechanicalParams * /*mparams*/ , Data& dOut, const Data& dIn) -{ - if( jacobian.rowSize() ) - { - auto dOutWa = sofa::helper::getWriteOnlyAccessor(dOut); - auto dInRa = sofa::helper::getReadAccessor(dIn); - jacobian.mult(dOutWa.wref(),dInRa.ref()); - } -} - -template -void SquareMapping::applyJT(const core::MechanicalParams * /*mparams*/ , Data& dIn, const Data& dOut) -{ - if( jacobian.rowSize() ) - { - auto dOutRa = sofa::helper::getReadAccessor(dOut); - auto dInWa = sofa::helper::getWriteOnlyAccessor(dIn); - jacobian.addMultTranspose(dInWa.wref(),dOutRa.ref()); - } + this->jacobian.compress(); } template -void SquareMapping::applyDJT(const core::MechanicalParams* mparams, core::MultiVecDerivId parentDfId, core::ConstMultiVecDerivId ) +void SquareMapping::matrixFreeApplyDJT( + const core::MechanicalParams* mparams, Real kFactor, + Data>& parentForce, + const Data>& parentDisplacement, + const Data>& childForce) { - const unsigned geometricStiffness = d_geometricStiffness.getValue().getSelectedId(); - if( !geometricStiffness ) return; + helper::WriteAccessor parentForceAccessor(parentForce); + helper::ReadAccessor parentDisplacementAccessor(parentDisplacement); + helper::ReadAccessor childForceAccessor(childForce); - helper::WriteAccessor > parentForce (*parentDfId[this->fromModel.get()].write()); - helper::ReadAccessor > parentDisplacement (*mparams->readDx(this->fromModel.get())); // parent displacement - SReal kfactor = mparams->kFactor(); - helper::ReadAccessor > childForce (*mparams->readF(this->toModel.get())); + const size_t size = parentDisplacementAccessor.size(); + kFactor *= 2.0; - if(d_useGeometricStiffnessMatrix.getValue() && K.compressedMatrix.nonZeros() ) + for(unsigned i=0; i -void SquareMapping::applyJT(const core::ConstraintParams* cparams, Data& out, const Data& in) -{ - SOFA_UNUSED(cparams); - auto childMatRa = sofa::helper::getReadAccessor(in); - auto parentMatWa = sofa::helper::getWriteAccessor(out); - addMultTransposeEigen(parentMatWa.wref(), jacobian.compressedMatrix, childMatRa.ref()); -} - - -template -const sofa::linearalgebra::BaseMatrix* SquareMapping::getJ() -{ - return &jacobian; } -template -const type::vector* SquareMapping::getJs() -{ - return &baseMatrices; -} - - - template void SquareMapping::updateK(const core::MechanicalParams *mparams, core::ConstMultiVecDerivId childForceId ) { 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(); + if( !geometricStiffness ) { this->K.resize(0,0); return; } - helper::ReadAccessor > childForce( *childForceId[this->toModel.get()].read() ); + helper::ReadAccessor childForce( *childForceId[this->toModel.get()].read() ); unsigned int size = this->fromModel->getSize(); - K.resizeBlocks(size,size); - K.reserve( size ); + this->K.resizeBlocks(size,size); + this->K.reserve( size ); for( size_t i=0 ; iK.beginRow(i); + this->K.insertBack( i, i, 2*childForce[i][0] ); } - K.compress(); -} - -template -const linearalgebra::BaseMatrix* SquareMapping::getK() -{ - return &K; + this->K.compress(); } template void SquareMapping::buildGeometricStiffnessMatrix( sofa::core::GeometricStiffnessMatrix* matrices) { - const unsigned geometricStiffness = d_geometricStiffness.getValue().getSelectedId(); + const unsigned geometricStiffness = this->d_geometricStiffness.getValue().getSelectedId(); if( !geometricStiffness ) { return; diff --git a/Sofa/Component/Mapping/NonLinear/tests/SquareMapping_test.cpp b/Sofa/Component/Mapping/NonLinear/tests/SquareMapping_test.cpp index ab49b58673c..66479a8c324 100644 --- a/Sofa/Component/Mapping/NonLinear/tests/SquareMapping_test.cpp +++ b/Sofa/Component/Mapping/NonLinear/tests/SquareMapping_test.cpp @@ -84,15 +84,5 @@ TYPED_TEST( SquareMappingTest , test ) ASSERT_TRUE(this->test()); } -TYPED_TEST( SquareMappingTest , testNoGeometricStiffnessMatrix ) -{ - TypeParam* map = static_cast( this->mapping ); - map->d_useGeometricStiffnessMatrix.setValue(false); - - ASSERT_TRUE(this->test()); -} - - - } // namespace } // namespace sofa From c3d12383cbb221a0b626ed2dff934a82bf5859ad Mon Sep 17 00:00:00 2001 From: Alex Bilger Date: Wed, 18 Sep 2024 16:29:26 +0200 Subject: [PATCH 06/16] delegate on updateK --- .../component/mapping/nonlinear/AreaMapping.h | 7 +++++- .../mapping/nonlinear/AreaMapping.inl | 18 ++++--------- .../nonlinear/AssembledNonLinearMapping.h | 11 ++++++-- .../nonlinear/AssembledNonLinearMapping.inl | 19 ++++++++++++++ .../nonlinear/DistanceFromTargetMapping.h | 7 +++++- .../nonlinear/DistanceFromTargetMapping.inl | 18 ++++++------- .../mapping/nonlinear/DistanceMapping.h | 7 +++++- .../mapping/nonlinear/DistanceMapping.inl | 25 ++++++++----------- .../mapping/nonlinear/SquareDistanceMapping.h | 7 +++++- .../nonlinear/SquareDistanceMapping.inl | 21 ++++++++-------- .../mapping/nonlinear/SquareMapping.h | 7 +++++- .../mapping/nonlinear/SquareMapping.inl | 17 ++++++------- .../mapping/nonlinear/VolumeMapping.h | 7 +++++- .../mapping/nonlinear/VolumeMapping.inl | 22 ++++++---------- 14 files changed, 112 insertions(+), 81 deletions(-) diff --git a/Sofa/Component/Mapping/NonLinear/src/sofa/component/mapping/nonlinear/AreaMapping.h b/Sofa/Component/Mapping/NonLinear/src/sofa/component/mapping/nonlinear/AreaMapping.h index 24457964f22..958e65dab29 100644 --- a/Sofa/Component/Mapping/NonLinear/src/sofa/component/mapping/nonlinear/AreaMapping.h +++ b/Sofa/Component/Mapping/NonLinear/src/sofa/component/mapping/nonlinear/AreaMapping.h @@ -51,7 +51,6 @@ class AreaMapping : public AssembledNonLinearMapping void init() override; void apply(const core::MechanicalParams* mparams, DataVecCoord_t& out, const DataVecCoord_t& in) override; - void updateK( const core::MechanicalParams* mparams, core::ConstMultiVecDerivId childForceId) override; void buildGeometricStiffnessMatrix(sofa::core::GeometricStiffnessMatrix* matrices) override; protected: @@ -62,6 +61,12 @@ class AreaMapping : public AssembledNonLinearMapping const Data >& parentDisplacement, const Data >& childForce) override; + using typename Inherit1::SparseKMatrixEigen; + + void doUpdateK( + const core::MechanicalParams* mparams, const Data >& childForce, + SparseKMatrixEigen& matrix) override; + const VecCoord_t* m_vertices{nullptr}; using Inherit1::JacobianEntry; diff --git a/Sofa/Component/Mapping/NonLinear/src/sofa/component/mapping/nonlinear/AreaMapping.inl b/Sofa/Component/Mapping/NonLinear/src/sofa/component/mapping/nonlinear/AreaMapping.inl index bc3f2c741fe..80502828dac 100644 --- a/Sofa/Component/Mapping/NonLinear/src/sofa/component/mapping/nonlinear/AreaMapping.inl +++ b/Sofa/Component/Mapping/NonLinear/src/sofa/component/mapping/nonlinear/AreaMapping.inl @@ -238,24 +238,18 @@ void AreaMapping::matrixFreeApplyDJT( } template -void AreaMapping::updateK(const core::MechanicalParams* mparams, - core::ConstMultiVecDerivId childForceId) +void AreaMapping::doUpdateK(const core::MechanicalParams* mparams, + const Data>& childForce, SparseKMatrixEigen& matrix) { SOFA_UNUSED(mparams); const unsigned geometricStiffness = this->d_geometricStiffness.getValue().getSelectedId(); - if( !geometricStiffness ) { this->K.resize(0,0); return; } - helper::ReadAccessor > > childForce( *childForceId[this->toModel.get()].read() ); - - { - unsigned int kSize = this->fromModel->getSize(); - this->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& childForceTri = childForce[triangleId]; + const Deriv_t& childForceTri = childForceAccessor[triangleId]; if( childForceTri[0] < 0 || geometricStiffness==1 ) { @@ -274,13 +268,11 @@ void AreaMapping::updateK(const core::MechanicalParams* mparams, { for (unsigned int j = 0; j < 3; ++j) { - this->K.addBlock(triangle[i], triangle[j], d2Area_d2x[i][j] * childForceTri[0]); + matrix.addBlock(triangle[i], triangle[j], d2Area_d2x[i][j] * childForceTri[0]); } } } } - - this->K.compress(); } template diff --git a/Sofa/Component/Mapping/NonLinear/src/sofa/component/mapping/nonlinear/AssembledNonLinearMapping.h b/Sofa/Component/Mapping/NonLinear/src/sofa/component/mapping/nonlinear/AssembledNonLinearMapping.h index 86990facad1..2bf10f76f3d 100644 --- a/Sofa/Component/Mapping/NonLinear/src/sofa/component/mapping/nonlinear/AssembledNonLinearMapping.h +++ b/Sofa/Component/Mapping/NonLinear/src/sofa/component/mapping/nonlinear/AssembledNonLinearMapping.h @@ -55,19 +55,26 @@ class AssembledNonLinearMapping : public core::Mapping, public NonLin const linearalgebra::BaseMatrix* getK() final; const type::vector* getJs() override; + void updateK( const core::MechanicalParams* mparams, core::ConstMultiVecDerivId childForceId) final; + protected: + + using SparseKMatrixEigen = linearalgebra::EigenSparseMatrix; + virtual void matrixFreeApplyDJT(const core::MechanicalParams* mparams, Real kFactor, Data >& parentForce, const Data >& parentDisplacement, const Data >& childForce) = 0; - using SparseKMatrixEigen = linearalgebra::EigenSparseMatrix; - SparseMatrixEigen jacobian; ///< Jacobian of the mapping type::vector baseMatrices; ///< Jacobian of the mapping, in a vector SparseKMatrixEigen K; ///< Assembled geometric stiffness matrix + virtual void doUpdateK( + const core::MechanicalParams* mparams, const Data >& childForce, + SparseKMatrixEigen& matrix) = 0; + /** * @brief Represents an entry in the Jacobian matrix. * diff --git a/Sofa/Component/Mapping/NonLinear/src/sofa/component/mapping/nonlinear/AssembledNonLinearMapping.inl b/Sofa/Component/Mapping/NonLinear/src/sofa/component/mapping/nonlinear/AssembledNonLinearMapping.inl index d636d1ccad7..39d46bf4595 100644 --- a/Sofa/Component/Mapping/NonLinear/src/sofa/component/mapping/nonlinear/AssembledNonLinearMapping.inl +++ b/Sofa/Component/Mapping/NonLinear/src/sofa/component/mapping/nonlinear/AssembledNonLinearMapping.inl @@ -114,5 +114,24 @@ const type::vector* AssembledNonLinearMapping< return &baseMatrices; } +template +void AssembledNonLinearMapping:: +updateK(const core::MechanicalParams* mparams, core::ConstMultiVecDerivId childForceId) +{ + const unsigned geometricStiffness = this->d_geometricStiffness.getValue().getSelectedId(); + if( !geometricStiffness ) { this->K.resize(0,0); return; } + + const Data >& childForce = *childForceId[this->toModel.get()].read(); + + { + unsigned int kSize = this->fromModel->getSize(); + K.resizeBlocks(kSize, kSize); + } + + doUpdateK(mparams, childForce, K); + + K.compress(); +} + } diff --git a/Sofa/Component/Mapping/NonLinear/src/sofa/component/mapping/nonlinear/DistanceFromTargetMapping.h b/Sofa/Component/Mapping/NonLinear/src/sofa/component/mapping/nonlinear/DistanceFromTargetMapping.h index 2218a80f9f2..fae60d5097a 100644 --- a/Sofa/Component/Mapping/NonLinear/src/sofa/component/mapping/nonlinear/DistanceFromTargetMapping.h +++ b/Sofa/Component/Mapping/NonLinear/src/sofa/component/mapping/nonlinear/DistanceFromTargetMapping.h @@ -104,7 +104,6 @@ class DistanceFromTargetMapping : public AssembledNonLinearMapping& out, const DataVecCoord_t& in) override; - void updateK( const core::MechanicalParams* mparams, core::ConstMultiVecDerivId childForce ) override; void buildGeometricStiffnessMatrix(sofa::core::GeometricStiffnessMatrix* matrices) override; void draw(const core::visual::VisualParams* vparams) override; @@ -120,6 +119,12 @@ class DistanceFromTargetMapping : public AssembledNonLinearMapping >& parentDisplacement, const Data >& childForce) override; + using typename Inherit1::SparseKMatrixEigen; + + void doUpdateK( + const core::MechanicalParams* mparams, const Data >& childForce, + SparseKMatrixEigen& matrix) override; + type::vector directions; ///< Unit vectors in the directions of the lines type::vector< Real > invlengths; ///< inverse of current distances. Null represents the infinity (null distance) diff --git a/Sofa/Component/Mapping/NonLinear/src/sofa/component/mapping/nonlinear/DistanceFromTargetMapping.inl b/Sofa/Component/Mapping/NonLinear/src/sofa/component/mapping/nonlinear/DistanceFromTargetMapping.inl index 1480c7b8301..c3e25d10967 100644 --- a/Sofa/Component/Mapping/NonLinear/src/sofa/component/mapping/nonlinear/DistanceFromTargetMapping.inl +++ b/Sofa/Component/Mapping/NonLinear/src/sofa/component/mapping/nonlinear/DistanceFromTargetMapping.inl @@ -271,23 +271,22 @@ void DistanceFromTargetMapping::buildGeometricStiffnessMatrix( } template -void DistanceFromTargetMapping::updateK( const core::MechanicalParams* mparams, core::ConstMultiVecDerivId childForceId ) +void DistanceFromTargetMapping::doUpdateK( + const core::MechanicalParams* mparams, + const Data>& childForce, SparseKMatrixEigen& matrix) { SOFA_UNUSED(mparams); const unsigned geometricStiffness = this->d_geometricStiffness.getValue().getSelectedId(); - if( !geometricStiffness ) { K.resize(0,0); return; } - helper::ReadAccessor childForce( *childForceId[this->toModel.get()].read() ); + const helper::ReadAccessor childForceAccessor(childForce); const helper::ReadAccessor indices(d_indices); - helper::ReadAccessor in (*this->fromModel->read(core::ConstVecCoordId::position())); - K.resizeBlocks(in.size(),in.size()); - for(size_t i=0; i0) can lead to negative eigen values in geometric stiffness // this results in an undefinite implicit matrix that causes instabilities // if stabilized GS (geometricStiffness==2) -> keep only force in extension - if( childForce[i][0] < 0 || geometricStiffness==1 ) + if( childForceAccessor[i][0] < 0 || geometricStiffness==1 ) { size_t idx = indices[i]; @@ -299,12 +298,11 @@ void DistanceFromTargetMapping::updateK( const core::MechanicalParams b[j][k] = static_cast(1) * ( j==k ) - directions[i][j]*directions[i][k]; } } - b *= childForce[i][0] * invlengths[i]; // (I - uu^T)*f/l + b *= childForceAccessor[i][0] * invlengths[i]; // (I - uu^T)*f/l - K.addBlock(idx,idx,b); + matrix.addBlock(idx,idx,b); } } - K.compress(); } diff --git a/Sofa/Component/Mapping/NonLinear/src/sofa/component/mapping/nonlinear/DistanceMapping.h b/Sofa/Component/Mapping/NonLinear/src/sofa/component/mapping/nonlinear/DistanceMapping.h index caee530a3cb..9f0debd8b40 100644 --- a/Sofa/Component/Mapping/NonLinear/src/sofa/component/mapping/nonlinear/DistanceMapping.h +++ b/Sofa/Component/Mapping/NonLinear/src/sofa/component/mapping/nonlinear/DistanceMapping.h @@ -79,7 +79,6 @@ class DistanceMapping : public AssembledNonLinearMapping void init() override; void apply(const core::MechanicalParams *mparams, DataVecCoord_t& out, const DataVecCoord_t& in) override; - void updateK( const core::MechanicalParams* mparams, core::ConstMultiVecDerivId childForce ) override; void buildGeometricStiffnessMatrix(sofa::core::GeometricStiffnessMatrix* matrices) override; void computeBBox(const core::ExecParams* params, bool onlyVisible) override; @@ -93,6 +92,12 @@ class DistanceMapping : public AssembledNonLinearMapping const Data >& parentDisplacement, const Data >& childForce) override; + using typename Inherit1::SparseKMatrixEigen; + + void doUpdateK( + const core::MechanicalParams* mparams, const Data >& childForce, + SparseKMatrixEigen& matrix) override; + type::vector directions; ///< Unit vectors in the directions of the lines type::vector< Real > invlengths; ///< inverse of current distances. Null represents the infinity (null distance) diff --git a/Sofa/Component/Mapping/NonLinear/src/sofa/component/mapping/nonlinear/DistanceMapping.inl b/Sofa/Component/Mapping/NonLinear/src/sofa/component/mapping/nonlinear/DistanceMapping.inl index 10c682f833f..a2096e0e05e 100644 --- a/Sofa/Component/Mapping/NonLinear/src/sofa/component/mapping/nonlinear/DistanceMapping.inl +++ b/Sofa/Component/Mapping/NonLinear/src/sofa/component/mapping/nonlinear/DistanceMapping.inl @@ -211,26 +211,24 @@ void DistanceMapping::matrixFreeApplyDJT( } } - template -void DistanceMapping::updateK(const core::MechanicalParams *mparams, core::ConstMultiVecDerivId childForceId ) +void DistanceMapping::doUpdateK( + const core::MechanicalParams* mparams, + const Data>& childForce, SparseKMatrixEigen& matrix) { SOFA_UNUSED(mparams); const unsigned geometricStiffness = this->d_geometricStiffness.getValue().getSelectedId(); - if( !geometricStiffness ) { this->K.resize(0,0); return; } - const helper::ReadAccessor childForce( *childForceId[this->toModel.get()].read() ); + const helper::ReadAccessor childForceAccessor(childForce); const SeqEdges& links = l_topology->getEdges(); - unsigned int size = this->fromModel->getSize(); - this->K.resizeBlocks(size,size); - for(size_t i=0; i0) can lead to negative eigen values in geometric stiffness // this results in an undefinite implicit matrix that causes instabilities // if stabilized GS (geometricStiffness==2) -> keep only force in extension - if( childForce[i][0] < 0 || geometricStiffness==1 ) + if( childForceAccessor[i][0] < 0 || geometricStiffness==1 ) { sofa::type::Mat b; // = (I - uu^T) @@ -241,16 +239,15 @@ void DistanceMapping::updateK(const core::MechanicalParams *mparams, b[j][k] = static_cast(1) * ( j==k ) - directions[i][j]*directions[i][k]; } } - b *= childForce[i][0] * invlengths[i]; // (I - uu^T)*f/l + b *= childForceAccessor[i][0] * invlengths[i]; // (I - uu^T)*f/l // Note that 'links' is not sorted so the matrix can not be filled-up in order - this->K.addBlock(links[i][0],links[i][0],b); - this->K.addBlock(links[i][0],links[i][1],-b); - this->K.addBlock(links[i][1],links[i][0],-b); - this->K.addBlock(links[i][1],links[i][1],b); + matrix.addBlock(links[i][0],links[i][0],b); + matrix.addBlock(links[i][0],links[i][1],-b); + matrix.addBlock(links[i][1],links[i][0],-b); + matrix.addBlock(links[i][1],links[i][1],b); } } - this->K.compress(); } template diff --git a/Sofa/Component/Mapping/NonLinear/src/sofa/component/mapping/nonlinear/SquareDistanceMapping.h b/Sofa/Component/Mapping/NonLinear/src/sofa/component/mapping/nonlinear/SquareDistanceMapping.h index 491cdbcb699..2a55fe4edc3 100644 --- a/Sofa/Component/Mapping/NonLinear/src/sofa/component/mapping/nonlinear/SquareDistanceMapping.h +++ b/Sofa/Component/Mapping/NonLinear/src/sofa/component/mapping/nonlinear/SquareDistanceMapping.h @@ -76,7 +76,6 @@ class SquareDistanceMapping : public AssembledNonLinearMapping void init() override; void apply(const core::MechanicalParams *mparams, DataVecCoord_t& out, const DataVecCoord_t& in) override; - void updateK( const core::MechanicalParams* mparams, core::ConstMultiVecDerivId childForce ) override; void buildGeometricStiffnessMatrix(sofa::core::GeometricStiffnessMatrix* matrices) override; void draw(const core::visual::VisualParams* vparams) override; @@ -90,6 +89,12 @@ class SquareDistanceMapping : public AssembledNonLinearMapping const Data >& parentDisplacement, const Data >& childForce) override; + using typename Inherit1::SparseKMatrixEigen; + + void doUpdateK( + const core::MechanicalParams* mparams, const Data >& childForce, + SparseKMatrixEigen& matrix) override; + /// r=b-a only for position (eventual rotation, affine transform... remains null) void computeCoordPositionDifference( Direction& r, const Coord_t& a, const Coord_t& b ); }; diff --git a/Sofa/Component/Mapping/NonLinear/src/sofa/component/mapping/nonlinear/SquareDistanceMapping.inl b/Sofa/Component/Mapping/NonLinear/src/sofa/component/mapping/nonlinear/SquareDistanceMapping.inl index a88a9e02647..7891e0862ba 100644 --- a/Sofa/Component/Mapping/NonLinear/src/sofa/component/mapping/nonlinear/SquareDistanceMapping.inl +++ b/Sofa/Component/Mapping/NonLinear/src/sofa/component/mapping/nonlinear/SquareDistanceMapping.inl @@ -178,36 +178,35 @@ void SquareDistanceMapping::matrixFreeApplyDJT( } template -void SquareDistanceMapping::updateK(const core::MechanicalParams *mparams, core::ConstMultiVecDerivId childForceId ) +void SquareDistanceMapping::doUpdateK( + const core::MechanicalParams* mparams, + const Data>& childForce, SparseKMatrixEigen& matrix) { SOFA_UNUSED(mparams); const unsigned geometricStiffness = this->d_geometricStiffness.getValue().getSelectedId(); - if( !geometricStiffness ) { this->K.resize(0,0); return; } - helper::ReadAccessor childForce( *childForceId[this->toModel.get()].read() ); + const helper::ReadAccessor childForceAccessor(childForce); const SeqEdges& links = l_topology->getEdges(); unsigned int size = this->fromModel->getSize(); - this->K.resizeBlocks(size,size); for(size_t i=0; i0) can lead to negative eigen values in geometric stiffness // this results in an undefinite implicit matrix that causes instabilities // if stabilized GS (geometricStiffness==2) -> keep only force in extension - if( childForce[i][0] < 0 || geometricStiffness==1 ) + if( childForceAccessor[i][0] < 0 || geometricStiffness==1 ) { - SReal tmp = 2*childForce[i][0]; + SReal tmp = 2*childForceAccessor[i][0]; for(unsigned k=0; kK.add( links[i][0]*Nin+k, links[i][0]*Nin+k, tmp ); - this->K.add( links[i][0]*Nin+k, links[i][1]*Nin+k, -tmp ); - this->K.add( links[i][1]*Nin+k, links[i][1]*Nin+k, tmp ); - this->K.add( links[i][1]*Nin+k, links[i][0]*Nin+k, -tmp ); + matrix.add( links[i][0]*Nin+k, links[i][0]*Nin+k, tmp ); + matrix.add( links[i][0]*Nin+k, links[i][1]*Nin+k, -tmp ); + matrix.add( links[i][1]*Nin+k, links[i][1]*Nin+k, tmp ); + matrix.add( links[i][1]*Nin+k, links[i][0]*Nin+k, -tmp ); } } } - this->K.compress(); } template diff --git a/Sofa/Component/Mapping/NonLinear/src/sofa/component/mapping/nonlinear/SquareMapping.h b/Sofa/Component/Mapping/NonLinear/src/sofa/component/mapping/nonlinear/SquareMapping.h index 40e89b9a14a..58130c043f1 100644 --- a/Sofa/Component/Mapping/NonLinear/src/sofa/component/mapping/nonlinear/SquareMapping.h +++ b/Sofa/Component/Mapping/NonLinear/src/sofa/component/mapping/nonlinear/SquareMapping.h @@ -53,7 +53,6 @@ class SquareMapping : public AssembledNonLinearMapping void init() override; void apply(const core::MechanicalParams *mparams, DataVecCoord_t& out, const DataVecCoord_t& in) override; - void updateK( const core::MechanicalParams* mparams, core::ConstMultiVecDerivId childForce ) override; void buildGeometricStiffnessMatrix(sofa::core::GeometricStiffnessMatrix* matrices) override; protected: @@ -62,6 +61,12 @@ class SquareMapping : public AssembledNonLinearMapping Data >& parentForce, const Data >& parentDisplacement, const Data >& childForce) override; + + using typename Inherit1::SparseKMatrixEigen; + + void doUpdateK( + const core::MechanicalParams* mparams, const Data >& childForce, + SparseKMatrixEigen& matrix) override; }; diff --git a/Sofa/Component/Mapping/NonLinear/src/sofa/component/mapping/nonlinear/SquareMapping.inl b/Sofa/Component/Mapping/NonLinear/src/sofa/component/mapping/nonlinear/SquareMapping.inl index 9f69a3893cd..e5a768862ad 100644 --- a/Sofa/Component/Mapping/NonLinear/src/sofa/component/mapping/nonlinear/SquareMapping.inl +++ b/Sofa/Component/Mapping/NonLinear/src/sofa/component/mapping/nonlinear/SquareMapping.inl @@ -88,23 +88,20 @@ void SquareMapping::matrixFreeApplyDJT( } template -void SquareMapping::updateK(const core::MechanicalParams *mparams, core::ConstMultiVecDerivId childForceId ) +void SquareMapping::doUpdateK(const core::MechanicalParams* mparams, + const Data>& childForce, SparseKMatrixEigen& matrix) { SOFA_UNUSED(mparams); const unsigned geometricStiffness = this->d_geometricStiffness.getValue().getSelectedId(); - if( !geometricStiffness ) { this->K.resize(0,0); return; } - - helper::ReadAccessor childForce( *childForceId[this->toModel.get()].read() ); + const helper::ReadAccessor childForceAccessor(childForce); unsigned int size = this->fromModel->getSize(); - this->K.resizeBlocks(size,size); - this->K.reserve( size ); - for( size_t i=0 ; iK.beginRow(i); - this->K.insertBack( i, i, 2*childForce[i][0] ); + matrix.beginRow(i); + matrix.insertBack( i, i, 2*childForceAccessor[i][0] ); } - this->K.compress(); } template diff --git a/Sofa/Component/Mapping/NonLinear/src/sofa/component/mapping/nonlinear/VolumeMapping.h b/Sofa/Component/Mapping/NonLinear/src/sofa/component/mapping/nonlinear/VolumeMapping.h index e677bb70503..d930beb2251 100644 --- a/Sofa/Component/Mapping/NonLinear/src/sofa/component/mapping/nonlinear/VolumeMapping.h +++ b/Sofa/Component/Mapping/NonLinear/src/sofa/component/mapping/nonlinear/VolumeMapping.h @@ -51,7 +51,6 @@ class VolumeMapping : public AssembledNonLinearMapping void init() override; void apply(const core::MechanicalParams* mparams, DataVecCoord_t& out, const DataVecCoord_t& in) override; - void updateK( const core::MechanicalParams* mparams, core::ConstMultiVecDerivId childForceId) override; void buildGeometricStiffnessMatrix(sofa::core::GeometricStiffnessMatrix* matrices) override; protected: @@ -62,6 +61,12 @@ class VolumeMapping : public AssembledNonLinearMapping const Data >& parentDisplacement, const Data >& childForce) override; + using typename Inherit1::SparseKMatrixEigen; + + void doUpdateK( + const core::MechanicalParams* mparams, const Data >& childForce, + SparseKMatrixEigen& matrix) override; + const VecCoord_t* m_vertices{nullptr}; using Inherit1::JacobianEntry; diff --git a/Sofa/Component/Mapping/NonLinear/src/sofa/component/mapping/nonlinear/VolumeMapping.inl b/Sofa/Component/Mapping/NonLinear/src/sofa/component/mapping/nonlinear/VolumeMapping.inl index 94e6d4b266f..62e53757919 100644 --- a/Sofa/Component/Mapping/NonLinear/src/sofa/component/mapping/nonlinear/VolumeMapping.inl +++ b/Sofa/Component/Mapping/NonLinear/src/sofa/component/mapping/nonlinear/VolumeMapping.inl @@ -241,26 +241,20 @@ void VolumeMapping::matrixFreeApplyDJT( } template -void VolumeMapping::updateK(const core::MechanicalParams* mparams, - core::ConstMultiVecDerivId childForceId) +void VolumeMapping::doUpdateK(const core::MechanicalParams* mparams, + const Data>& childForce, SparseKMatrixEigen& matrix) { SOFA_UNUSED(mparams); const unsigned geometricStiffness = this->d_geometricStiffness.getValue().getSelectedId(); - if( !geometricStiffness ) { this->K.resize(0,0); return; } - helper::ReadAccessor > > childForce( *childForceId[this->toModel.get()].read() ); - - { - unsigned int kSize = this->fromModel->getSize(); - this->K.resizeBlocks(kSize, kSize); - } + const helper::ReadAccessor childForceAccessor(childForce); const auto& tetrahedra = l_topology->getTetrahedra(); for (unsigned int tetId = 0; tetId < tetrahedra.size(); ++tetId) { - const Deriv_t& childForceTri = childForce[tetId]; + const Deriv_t& childForceTetra = childForceAccessor[tetId]; - if( childForceTri[0] < 0 || geometricStiffness==1 ) + if( childForceTetra[0] < 0 || geometricStiffness==1 ) { const auto& tetra = tetrahedra[tetId]; @@ -278,14 +272,12 @@ void VolumeMapping::updateK(const core::MechanicalParams* mparams, { for (unsigned int j = i+1; j < 4; ++j) //diagonal terms are omitted because they are null { - this->K.addBlock(tetra[i], tetra[j], d2Volume_d2x(i, j) * childForceTri[0]); - this->K.addBlock(tetra[j], tetra[i], d2Volume_d2x(j, i) * childForceTri[0]); + matrix.addBlock(tetra[i], tetra[j], d2Volume_d2x(i, j) * childForceTetra[0]); + matrix.addBlock(tetra[j], tetra[i], d2Volume_d2x(j, i) * childForceTetra[0]); } } } } - - this->K.compress(); } template From eced036ccce7901ade7c93e7172c8ea5f106f4da Mon Sep 17 00:00:00 2001 From: Alex Bilger Date: Wed, 18 Sep 2024 16:31:04 +0200 Subject: [PATCH 07/16] K can be private --- .../component/mapping/nonlinear/AssembledNonLinearMapping.h | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/Sofa/Component/Mapping/NonLinear/src/sofa/component/mapping/nonlinear/AssembledNonLinearMapping.h b/Sofa/Component/Mapping/NonLinear/src/sofa/component/mapping/nonlinear/AssembledNonLinearMapping.h index 2bf10f76f3d..ac79241fc41 100644 --- a/Sofa/Component/Mapping/NonLinear/src/sofa/component/mapping/nonlinear/AssembledNonLinearMapping.h +++ b/Sofa/Component/Mapping/NonLinear/src/sofa/component/mapping/nonlinear/AssembledNonLinearMapping.h @@ -69,7 +69,6 @@ class AssembledNonLinearMapping : public core::Mapping, public NonLin SparseMatrixEigen jacobian; ///< Jacobian of the mapping type::vector baseMatrices; ///< Jacobian of the mapping, in a vector - SparseKMatrixEigen K; ///< Assembled geometric stiffness matrix virtual void doUpdateK( const core::MechanicalParams* mparams, const Data >& childForce, @@ -89,6 +88,10 @@ class AssembledNonLinearMapping : public core::Mapping, public NonLin typename In::CPos jacobianValue; bool operator<(const JacobianEntry& other) const { return vertexId < other.vertexId;} }; + +private: + + SparseKMatrixEigen K; ///< Assembled geometric stiffness matrix }; } From fc0ca98d28b950b99c4c4c04949e1c6648f46620 Mon Sep 17 00:00:00 2001 From: Alex Bilger Date: Wed, 18 Sep 2024 16:33:38 +0200 Subject: [PATCH 08/16] factorize baseMatrix --- .../sofa/component/mapping/nonlinear/AreaMapping.inl | 3 --- .../mapping/nonlinear/AssembledNonLinearMapping.inl | 3 +++ .../mapping/nonlinear/DistanceFromTargetMapping.inl | 10 +++------- .../component/mapping/nonlinear/DistanceMapping.inl | 3 --- .../mapping/nonlinear/SquareDistanceMapping.inl | 3 --- .../sofa/component/mapping/nonlinear/SquareMapping.h | 2 -- .../sofa/component/mapping/nonlinear/SquareMapping.inl | 9 --------- .../sofa/component/mapping/nonlinear/VolumeMapping.inl | 3 --- 8 files changed, 6 insertions(+), 30 deletions(-) diff --git a/Sofa/Component/Mapping/NonLinear/src/sofa/component/mapping/nonlinear/AreaMapping.inl b/Sofa/Component/Mapping/NonLinear/src/sofa/component/mapping/nonlinear/AreaMapping.inl index 80502828dac..1d85b69b91d 100644 --- a/Sofa/Component/Mapping/NonLinear/src/sofa/component/mapping/nonlinear/AreaMapping.inl +++ b/Sofa/Component/Mapping/NonLinear/src/sofa/component/mapping/nonlinear/AreaMapping.inl @@ -125,9 +125,6 @@ void AreaMapping::init() this->getToModel()->resize( nbTriangles ); this->jacobian.resizeBlocks(nbTriangles, pos.size()); - this->baseMatrices.resize( 1 ); - this->baseMatrices[0] = &this->jacobian; - Inherit1::init(); if (this->d_componentState.getValue() != sofa::core::objectmodel::ComponentState::Invalid) diff --git a/Sofa/Component/Mapping/NonLinear/src/sofa/component/mapping/nonlinear/AssembledNonLinearMapping.inl b/Sofa/Component/Mapping/NonLinear/src/sofa/component/mapping/nonlinear/AssembledNonLinearMapping.inl index 39d46bf4595..47417e63d13 100644 --- a/Sofa/Component/Mapping/NonLinear/src/sofa/component/mapping/nonlinear/AssembledNonLinearMapping.inl +++ b/Sofa/Component/Mapping/NonLinear/src/sofa/component/mapping/nonlinear/AssembledNonLinearMapping.inl @@ -32,6 +32,9 @@ void AssembledNonLinearMapping:: init() { core::Mapping::init(); + + this->baseMatrices.resize( 1 ); + this->baseMatrices[0] = &this->jacobian; } template diff --git a/Sofa/Component/Mapping/NonLinear/src/sofa/component/mapping/nonlinear/DistanceFromTargetMapping.inl b/Sofa/Component/Mapping/NonLinear/src/sofa/component/mapping/nonlinear/DistanceFromTargetMapping.inl index c3e25d10967..31cf683e576 100644 --- a/Sofa/Component/Mapping/NonLinear/src/sofa/component/mapping/nonlinear/DistanceFromTargetMapping.inl +++ b/Sofa/Component/Mapping/NonLinear/src/sofa/component/mapping/nonlinear/DistanceFromTargetMapping.inl @@ -118,10 +118,6 @@ void DistanceFromTargetMapping::init() this->getToModel()->resize(d_indices.getValue().size() ); - - baseMatrices.resize( 1 ); - baseMatrices[0] = &jacobian; - Inherit1::init(); // applies the mapping, so after the Data init } @@ -170,15 +166,15 @@ void DistanceFromTargetMapping::apply(const core::MechanicalParams * for(unsigned j=0; jjacobian.beginRow(i*Nout+j); for(unsigned k=0; kjacobian.insertBack( i*Nout+j, indices[i]*Nin+k, gap[k] ); } } } - jacobian.compress(); + this->jacobian.compress(); } template diff --git a/Sofa/Component/Mapping/NonLinear/src/sofa/component/mapping/nonlinear/DistanceMapping.inl b/Sofa/Component/Mapping/NonLinear/src/sofa/component/mapping/nonlinear/DistanceMapping.inl index a2096e0e05e..8e2a75fcdfe 100644 --- a/Sofa/Component/Mapping/NonLinear/src/sofa/component/mapping/nonlinear/DistanceMapping.inl +++ b/Sofa/Component/Mapping/NonLinear/src/sofa/component/mapping/nonlinear/DistanceMapping.inl @@ -91,9 +91,6 @@ void DistanceMapping::init() } } - this->baseMatrices.resize( 1 ); - this->baseMatrices[0] = &this->jacobian; - this->Inherit1::init(); // applies the mapping, so after the Data init } diff --git a/Sofa/Component/Mapping/NonLinear/src/sofa/component/mapping/nonlinear/SquareDistanceMapping.inl b/Sofa/Component/Mapping/NonLinear/src/sofa/component/mapping/nonlinear/SquareDistanceMapping.inl index 7891e0862ba..db6447dd0cb 100644 --- a/Sofa/Component/Mapping/NonLinear/src/sofa/component/mapping/nonlinear/SquareDistanceMapping.inl +++ b/Sofa/Component/Mapping/NonLinear/src/sofa/component/mapping/nonlinear/SquareDistanceMapping.inl @@ -69,9 +69,6 @@ void SquareDistanceMapping::init() this->getToModel()->resize( links.size() ); - this->baseMatrices.resize( 1 ); - this->baseMatrices[0] = &this->jacobian; - this->Inherit1::init(); } diff --git a/Sofa/Component/Mapping/NonLinear/src/sofa/component/mapping/nonlinear/SquareMapping.h b/Sofa/Component/Mapping/NonLinear/src/sofa/component/mapping/nonlinear/SquareMapping.h index 58130c043f1..86b7c51d867 100644 --- a/Sofa/Component/Mapping/NonLinear/src/sofa/component/mapping/nonlinear/SquareMapping.h +++ b/Sofa/Component/Mapping/NonLinear/src/sofa/component/mapping/nonlinear/SquareMapping.h @@ -50,8 +50,6 @@ class SquareMapping : public AssembledNonLinearMapping static constexpr auto Nin = In::deriv_total_size; - void init() override; - void apply(const core::MechanicalParams *mparams, DataVecCoord_t& out, const DataVecCoord_t& in) override; void buildGeometricStiffnessMatrix(sofa::core::GeometricStiffnessMatrix* matrices) override; diff --git a/Sofa/Component/Mapping/NonLinear/src/sofa/component/mapping/nonlinear/SquareMapping.inl b/Sofa/Component/Mapping/NonLinear/src/sofa/component/mapping/nonlinear/SquareMapping.inl index e5a768862ad..46486ce1493 100644 --- a/Sofa/Component/Mapping/NonLinear/src/sofa/component/mapping/nonlinear/SquareMapping.inl +++ b/Sofa/Component/Mapping/NonLinear/src/sofa/component/mapping/nonlinear/SquareMapping.inl @@ -33,15 +33,6 @@ namespace sofa::component::mapping::nonlinear { -template -void SquareMapping::init() -{ - this->baseMatrices.resize( 1 ); - this->baseMatrices[0] = &this->jacobian; - - Inherit1::init(); -} - template void SquareMapping::apply(const core::MechanicalParams* mparams, DataVecCoord_t& dOut, const DataVecCoord_t& dIn) diff --git a/Sofa/Component/Mapping/NonLinear/src/sofa/component/mapping/nonlinear/VolumeMapping.inl b/Sofa/Component/Mapping/NonLinear/src/sofa/component/mapping/nonlinear/VolumeMapping.inl index 62e53757919..c3ce47d3479 100644 --- a/Sofa/Component/Mapping/NonLinear/src/sofa/component/mapping/nonlinear/VolumeMapping.inl +++ b/Sofa/Component/Mapping/NonLinear/src/sofa/component/mapping/nonlinear/VolumeMapping.inl @@ -104,9 +104,6 @@ void VolumeMapping::init() this->getToModel()->resize( nbTetrahedra ); this->jacobian.resizeBlocks(nbTetrahedra, pos.size()); - this->baseMatrices.resize( 1 ); - this->baseMatrices[0] = &this->jacobian; - Inherit1::init(); if (this->d_componentState.getValue() != sofa::core::objectmodel::ComponentState::Invalid) From 41c7649507e221e8b2122dc017f2e161569f5ccb Mon Sep 17 00:00:00 2001 From: Alex Bilger Date: Wed, 18 Sep 2024 16:34:26 +0200 Subject: [PATCH 09/16] baseMatrices can be private --- .../component/mapping/nonlinear/AssembledNonLinearMapping.h | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Sofa/Component/Mapping/NonLinear/src/sofa/component/mapping/nonlinear/AssembledNonLinearMapping.h b/Sofa/Component/Mapping/NonLinear/src/sofa/component/mapping/nonlinear/AssembledNonLinearMapping.h index ac79241fc41..65e6333da1c 100644 --- a/Sofa/Component/Mapping/NonLinear/src/sofa/component/mapping/nonlinear/AssembledNonLinearMapping.h +++ b/Sofa/Component/Mapping/NonLinear/src/sofa/component/mapping/nonlinear/AssembledNonLinearMapping.h @@ -68,7 +68,6 @@ class AssembledNonLinearMapping : public core::Mapping, public NonLin const Data >& childForce) = 0; SparseMatrixEigen jacobian; ///< Jacobian of the mapping - type::vector baseMatrices; ///< Jacobian of the mapping, in a vector virtual void doUpdateK( const core::MechanicalParams* mparams, const Data >& childForce, @@ -92,6 +91,7 @@ class AssembledNonLinearMapping : public core::Mapping, public NonLin private: SparseKMatrixEigen K; ///< Assembled geometric stiffness matrix + type::vector baseMatrices; ///< Jacobian of the mapping, in a vector }; } From 2a95663170f6bb81aa5d208aef67c1f7c79b2464 Mon Sep 17 00:00:00 2001 From: Alex Bilger Date: Wed, 18 Sep 2024 16:36:15 +0200 Subject: [PATCH 10/16] rename base class --- .../component/mapping/nonlinear/AreaMapping.h | 4 ++-- .../nonlinear/AssembledNonLinearMapping.h | 4 ++-- .../nonlinear/AssembledNonLinearMapping.inl | 16 ++++++++-------- .../nonlinear/DistanceFromTargetMapping.h | 4 ++-- .../mapping/nonlinear/DistanceMapping.h | 4 ++-- .../mapping/nonlinear/SquareDistanceMapping.h | 4 ++-- .../component/mapping/nonlinear/SquareMapping.h | 4 ++-- .../component/mapping/nonlinear/VolumeMapping.h | 4 ++-- 8 files changed, 22 insertions(+), 22 deletions(-) diff --git a/Sofa/Component/Mapping/NonLinear/src/sofa/component/mapping/nonlinear/AreaMapping.h b/Sofa/Component/Mapping/NonLinear/src/sofa/component/mapping/nonlinear/AreaMapping.h index 958e65dab29..5ee64ed8d88 100644 --- a/Sofa/Component/Mapping/NonLinear/src/sofa/component/mapping/nonlinear/AreaMapping.h +++ b/Sofa/Component/Mapping/NonLinear/src/sofa/component/mapping/nonlinear/AreaMapping.h @@ -32,10 +32,10 @@ namespace sofa::component::mapping::nonlinear { template -class AreaMapping : public AssembledNonLinearMapping +class AreaMapping : public BaseNonLinearMapping { public: - SOFA_CLASS(SOFA_TEMPLATE2(AreaMapping,TIn,TOut), SOFA_TEMPLATE3(AssembledNonLinearMapping,TIn,TOut, true)); + SOFA_CLASS(SOFA_TEMPLATE2(AreaMapping,TIn,TOut), SOFA_TEMPLATE3(BaseNonLinearMapping,TIn,TOut, true)); using In = TIn; using Out = TOut; diff --git a/Sofa/Component/Mapping/NonLinear/src/sofa/component/mapping/nonlinear/AssembledNonLinearMapping.h b/Sofa/Component/Mapping/NonLinear/src/sofa/component/mapping/nonlinear/AssembledNonLinearMapping.h index 65e6333da1c..3df90b304ce 100644 --- a/Sofa/Component/Mapping/NonLinear/src/sofa/component/mapping/nonlinear/AssembledNonLinearMapping.h +++ b/Sofa/Component/Mapping/NonLinear/src/sofa/component/mapping/nonlinear/AssembledNonLinearMapping.h @@ -30,11 +30,11 @@ namespace sofa::component::mapping::nonlinear { template -class AssembledNonLinearMapping : public core::Mapping, public NonLinearMappingData +class BaseNonLinearMapping : public core::Mapping, public NonLinearMappingData { public: SOFA_CLASS( - SOFA_TEMPLATE3(AssembledNonLinearMapping,TIn,TOut, HasStabilizedGeometricStiffness), + SOFA_TEMPLATE3(BaseNonLinearMapping,TIn,TOut, HasStabilizedGeometricStiffness), SOFA_TEMPLATE2(core::Mapping,TIn,TOut)); using In = TIn; diff --git a/Sofa/Component/Mapping/NonLinear/src/sofa/component/mapping/nonlinear/AssembledNonLinearMapping.inl b/Sofa/Component/Mapping/NonLinear/src/sofa/component/mapping/nonlinear/AssembledNonLinearMapping.inl index 47417e63d13..c6cc5ba0c9e 100644 --- a/Sofa/Component/Mapping/NonLinear/src/sofa/component/mapping/nonlinear/AssembledNonLinearMapping.inl +++ b/Sofa/Component/Mapping/NonLinear/src/sofa/component/mapping/nonlinear/AssembledNonLinearMapping.inl @@ -28,7 +28,7 @@ namespace sofa::component::mapping::nonlinear { template -void AssembledNonLinearMapping:: +void BaseNonLinearMapping:: init() { core::Mapping::init(); @@ -38,7 +38,7 @@ init() } template -void AssembledNonLinearMapping::applyJ( +void BaseNonLinearMapping::applyJ( const core::MechanicalParams* mparams, DataVecDeriv_t& out, const DataVecDeriv_t& in) { @@ -51,7 +51,7 @@ void AssembledNonLinearMapping::appl } template -void AssembledNonLinearMapping::applyJT( +void BaseNonLinearMapping::applyJT( const core::MechanicalParams* mparams, DataVecDeriv_t& out, const DataVecDeriv_t& in) { @@ -64,7 +64,7 @@ void AssembledNonLinearMapping::appl } template -void AssembledNonLinearMapping::applyJT( +void BaseNonLinearMapping::applyJT( const core::ConstraintParams* cparams, DataMatrixDeriv_t& out, const DataMatrixDeriv_t& in) { @@ -75,7 +75,7 @@ void AssembledNonLinearMapping::appl } template -void AssembledNonLinearMapping::applyDJT( +void BaseNonLinearMapping::applyDJT( const core::MechanicalParams* mparams, core::MultiVecDerivId parentForceId, core::ConstMultiVecDerivId childForceId) { @@ -106,19 +106,19 @@ void AssembledNonLinearMapping::appl } template -const linearalgebra::BaseMatrix* AssembledNonLinearMapping::getK() +const linearalgebra::BaseMatrix* BaseNonLinearMapping::getK() { return &K; } template -const type::vector* AssembledNonLinearMapping::getJs() +const type::vector* BaseNonLinearMapping::getJs() { return &baseMatrices; } template -void AssembledNonLinearMapping:: +void BaseNonLinearMapping:: updateK(const core::MechanicalParams* mparams, core::ConstMultiVecDerivId childForceId) { const unsigned geometricStiffness = this->d_geometricStiffness.getValue().getSelectedId(); diff --git a/Sofa/Component/Mapping/NonLinear/src/sofa/component/mapping/nonlinear/DistanceFromTargetMapping.h b/Sofa/Component/Mapping/NonLinear/src/sofa/component/mapping/nonlinear/DistanceFromTargetMapping.h index fae60d5097a..8d7d3c01676 100644 --- a/Sofa/Component/Mapping/NonLinear/src/sofa/component/mapping/nonlinear/DistanceFromTargetMapping.h +++ b/Sofa/Component/Mapping/NonLinear/src/sofa/component/mapping/nonlinear/DistanceFromTargetMapping.h @@ -60,10 +60,10 @@ struct BaseDistanceFromTargetMapping * @author Francois Faure */ template -class DistanceFromTargetMapping : public AssembledNonLinearMapping, public BaseDistanceFromTargetMapping +class DistanceFromTargetMapping : public BaseNonLinearMapping, public BaseDistanceFromTargetMapping { public: - SOFA_CLASS(SOFA_TEMPLATE2(DistanceFromTargetMapping,TIn,TOut), SOFA_TEMPLATE3(AssembledNonLinearMapping,TIn,TOut, true)); + SOFA_CLASS(SOFA_TEMPLATE2(DistanceFromTargetMapping,TIn,TOut), SOFA_TEMPLATE3(BaseNonLinearMapping,TIn,TOut, true)); using In = TIn; using Out = TOut; diff --git a/Sofa/Component/Mapping/NonLinear/src/sofa/component/mapping/nonlinear/DistanceMapping.h b/Sofa/Component/Mapping/NonLinear/src/sofa/component/mapping/nonlinear/DistanceMapping.h index 9f0debd8b40..aa140874632 100644 --- a/Sofa/Component/Mapping/NonLinear/src/sofa/component/mapping/nonlinear/DistanceMapping.h +++ b/Sofa/Component/Mapping/NonLinear/src/sofa/component/mapping/nonlinear/DistanceMapping.h @@ -46,10 +46,10 @@ namespace sofa::component::mapping::nonlinear * @author Francois Faure */ template -class DistanceMapping : public AssembledNonLinearMapping +class DistanceMapping : public BaseNonLinearMapping { public: - SOFA_CLASS(SOFA_TEMPLATE2(DistanceMapping,TIn,TOut), SOFA_TEMPLATE3(AssembledNonLinearMapping,TIn,TOut,true)); + SOFA_CLASS(SOFA_TEMPLATE2(DistanceMapping,TIn,TOut), SOFA_TEMPLATE3(BaseNonLinearMapping,TIn,TOut,true)); typedef TIn In; typedef TOut Out; diff --git a/Sofa/Component/Mapping/NonLinear/src/sofa/component/mapping/nonlinear/SquareDistanceMapping.h b/Sofa/Component/Mapping/NonLinear/src/sofa/component/mapping/nonlinear/SquareDistanceMapping.h index 2a55fe4edc3..423e4dea316 100644 --- a/Sofa/Component/Mapping/NonLinear/src/sofa/component/mapping/nonlinear/SquareDistanceMapping.h +++ b/Sofa/Component/Mapping/NonLinear/src/sofa/component/mapping/nonlinear/SquareDistanceMapping.h @@ -52,10 +52,10 @@ namespace sofa::component::mapping::nonlinear // If the rest lengths are not defined, they are set using the initial values. // If computeDistance is set to true, the rest lengths are set to 0. template -class SquareDistanceMapping : public AssembledNonLinearMapping +class SquareDistanceMapping : public BaseNonLinearMapping { public: - SOFA_CLASS(SOFA_TEMPLATE2(SquareDistanceMapping,TIn,TOut), SOFA_TEMPLATE3(AssembledNonLinearMapping,TIn,TOut,true)); + SOFA_CLASS(SOFA_TEMPLATE2(SquareDistanceMapping,TIn,TOut), SOFA_TEMPLATE3(BaseNonLinearMapping,TIn,TOut,true)); using In = TIn; using Out = TOut; diff --git a/Sofa/Component/Mapping/NonLinear/src/sofa/component/mapping/nonlinear/SquareMapping.h b/Sofa/Component/Mapping/NonLinear/src/sofa/component/mapping/nonlinear/SquareMapping.h index 86b7c51d867..5c6dc7f701f 100644 --- a/Sofa/Component/Mapping/NonLinear/src/sofa/component/mapping/nonlinear/SquareMapping.h +++ b/Sofa/Component/Mapping/NonLinear/src/sofa/component/mapping/nonlinear/SquareMapping.h @@ -38,10 +38,10 @@ namespace sofa::component::mapping::nonlinear */ template -class SquareMapping : public AssembledNonLinearMapping +class SquareMapping : public BaseNonLinearMapping { public: - SOFA_CLASS(SOFA_TEMPLATE2(SquareMapping,TIn,TOut), SOFA_TEMPLATE3(AssembledNonLinearMapping,TIn,TOut,false)); + SOFA_CLASS(SOFA_TEMPLATE2(SquareMapping,TIn,TOut), SOFA_TEMPLATE3(BaseNonLinearMapping,TIn,TOut,false)); using In = TIn; using Out = TOut; diff --git a/Sofa/Component/Mapping/NonLinear/src/sofa/component/mapping/nonlinear/VolumeMapping.h b/Sofa/Component/Mapping/NonLinear/src/sofa/component/mapping/nonlinear/VolumeMapping.h index d930beb2251..debd8766eaa 100644 --- a/Sofa/Component/Mapping/NonLinear/src/sofa/component/mapping/nonlinear/VolumeMapping.h +++ b/Sofa/Component/Mapping/NonLinear/src/sofa/component/mapping/nonlinear/VolumeMapping.h @@ -32,10 +32,10 @@ namespace sofa::component::mapping::nonlinear { template -class VolumeMapping : public AssembledNonLinearMapping +class VolumeMapping : public BaseNonLinearMapping { public: - SOFA_CLASS(SOFA_TEMPLATE2(VolumeMapping,TIn,TOut), SOFA_TEMPLATE3(AssembledNonLinearMapping,TIn,TOut, true)); + SOFA_CLASS(SOFA_TEMPLATE2(VolumeMapping,TIn,TOut), SOFA_TEMPLATE3(BaseNonLinearMapping,TIn,TOut, true)); using In = TIn; using Out = TOut; From 5cf352d0773f1f47bf20a1b5798eabfbc33ca43f Mon Sep 17 00:00:00 2001 From: Alex Bilger Date: Wed, 18 Sep 2024 16:39:31 +0200 Subject: [PATCH 11/16] rename files --- Sofa/Component/Mapping/NonLinear/CMakeLists.txt | 4 ++-- .../src/sofa/component/mapping/nonlinear/AreaMapping.h | 2 +- .../src/sofa/component/mapping/nonlinear/AreaMapping.inl | 2 +- .../{AssembledNonLinearMapping.h => BaseNonLinearMapping.h} | 0 ...AssembledNonLinearMapping.inl => BaseNonLinearMapping.inl} | 2 +- .../component/mapping/nonlinear/DistanceFromTargetMapping.h | 2 +- .../component/mapping/nonlinear/DistanceFromTargetMapping.inl | 2 +- .../src/sofa/component/mapping/nonlinear/DistanceMapping.h | 2 +- .../src/sofa/component/mapping/nonlinear/DistanceMapping.inl | 2 +- .../sofa/component/mapping/nonlinear/SquareDistanceMapping.h | 2 +- .../component/mapping/nonlinear/SquareDistanceMapping.inl | 2 +- .../src/sofa/component/mapping/nonlinear/SquareMapping.h | 2 +- .../src/sofa/component/mapping/nonlinear/SquareMapping.inl | 2 +- .../src/sofa/component/mapping/nonlinear/VolumeMapping.h | 2 +- .../src/sofa/component/mapping/nonlinear/VolumeMapping.inl | 2 +- 15 files changed, 15 insertions(+), 15 deletions(-) rename Sofa/Component/Mapping/NonLinear/src/sofa/component/mapping/nonlinear/{AssembledNonLinearMapping.h => BaseNonLinearMapping.h} (100%) rename Sofa/Component/Mapping/NonLinear/src/sofa/component/mapping/nonlinear/{AssembledNonLinearMapping.inl => BaseNonLinearMapping.inl} (98%) diff --git a/Sofa/Component/Mapping/NonLinear/CMakeLists.txt b/Sofa/Component/Mapping/NonLinear/CMakeLists.txt index 03545b647ce..d97638d9c41 100644 --- a/Sofa/Component/Mapping/NonLinear/CMakeLists.txt +++ b/Sofa/Component/Mapping/NonLinear/CMakeLists.txt @@ -8,8 +8,8 @@ set(HEADER_FILES ${SOFACOMPONENTMAPPINGNONLINEAR_SOURCE_DIR}/init.h ${SOFACOMPONENTMAPPINGNONLINEAR_SOURCE_DIR}/AreaMapping.h ${SOFACOMPONENTMAPPINGNONLINEAR_SOURCE_DIR}/AreaMapping.inl - ${SOFACOMPONENTMAPPINGNONLINEAR_SOURCE_DIR}/AssembledNonLinearMapping.h - ${SOFACOMPONENTMAPPINGNONLINEAR_SOURCE_DIR}/AssembledNonLinearMapping.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 diff --git a/Sofa/Component/Mapping/NonLinear/src/sofa/component/mapping/nonlinear/AreaMapping.h b/Sofa/Component/Mapping/NonLinear/src/sofa/component/mapping/nonlinear/AreaMapping.h index 5ee64ed8d88..2392f310364 100644 --- a/Sofa/Component/Mapping/NonLinear/src/sofa/component/mapping/nonlinear/AreaMapping.h +++ b/Sofa/Component/Mapping/NonLinear/src/sofa/component/mapping/nonlinear/AreaMapping.h @@ -22,7 +22,7 @@ #pragma once #include -#include +#include #include #include #include diff --git a/Sofa/Component/Mapping/NonLinear/src/sofa/component/mapping/nonlinear/AreaMapping.inl b/Sofa/Component/Mapping/NonLinear/src/sofa/component/mapping/nonlinear/AreaMapping.inl index 1d85b69b91d..0a343e2fc2e 100644 --- a/Sofa/Component/Mapping/NonLinear/src/sofa/component/mapping/nonlinear/AreaMapping.inl +++ b/Sofa/Component/Mapping/NonLinear/src/sofa/component/mapping/nonlinear/AreaMapping.inl @@ -21,7 +21,7 @@ ******************************************************************************/ #pragma once #include -#include +#include #include namespace sofa::component::mapping::nonlinear diff --git a/Sofa/Component/Mapping/NonLinear/src/sofa/component/mapping/nonlinear/AssembledNonLinearMapping.h b/Sofa/Component/Mapping/NonLinear/src/sofa/component/mapping/nonlinear/BaseNonLinearMapping.h similarity index 100% rename from Sofa/Component/Mapping/NonLinear/src/sofa/component/mapping/nonlinear/AssembledNonLinearMapping.h rename to Sofa/Component/Mapping/NonLinear/src/sofa/component/mapping/nonlinear/BaseNonLinearMapping.h diff --git a/Sofa/Component/Mapping/NonLinear/src/sofa/component/mapping/nonlinear/AssembledNonLinearMapping.inl b/Sofa/Component/Mapping/NonLinear/src/sofa/component/mapping/nonlinear/BaseNonLinearMapping.inl similarity index 98% rename from Sofa/Component/Mapping/NonLinear/src/sofa/component/mapping/nonlinear/AssembledNonLinearMapping.inl rename to Sofa/Component/Mapping/NonLinear/src/sofa/component/mapping/nonlinear/BaseNonLinearMapping.inl index c6cc5ba0c9e..fec4f24a235 100644 --- a/Sofa/Component/Mapping/NonLinear/src/sofa/component/mapping/nonlinear/AssembledNonLinearMapping.inl +++ b/Sofa/Component/Mapping/NonLinear/src/sofa/component/mapping/nonlinear/BaseNonLinearMapping.inl @@ -21,7 +21,7 @@ ******************************************************************************/ #pragma once -#include +#include #include namespace sofa::component::mapping::nonlinear diff --git a/Sofa/Component/Mapping/NonLinear/src/sofa/component/mapping/nonlinear/DistanceFromTargetMapping.h b/Sofa/Component/Mapping/NonLinear/src/sofa/component/mapping/nonlinear/DistanceFromTargetMapping.h index 8d7d3c01676..9c226f13918 100644 --- a/Sofa/Component/Mapping/NonLinear/src/sofa/component/mapping/nonlinear/DistanceFromTargetMapping.h +++ b/Sofa/Component/Mapping/NonLinear/src/sofa/component/mapping/nonlinear/DistanceFromTargetMapping.h @@ -22,7 +22,7 @@ #pragma once #include -#include +#include #include #include #include diff --git a/Sofa/Component/Mapping/NonLinear/src/sofa/component/mapping/nonlinear/DistanceFromTargetMapping.inl b/Sofa/Component/Mapping/NonLinear/src/sofa/component/mapping/nonlinear/DistanceFromTargetMapping.inl index 31cf683e576..a5fb0255dc6 100644 --- a/Sofa/Component/Mapping/NonLinear/src/sofa/component/mapping/nonlinear/DistanceFromTargetMapping.inl +++ b/Sofa/Component/Mapping/NonLinear/src/sofa/component/mapping/nonlinear/DistanceFromTargetMapping.inl @@ -22,7 +22,7 @@ #pragma once #include -#include +#include #include #include #include diff --git a/Sofa/Component/Mapping/NonLinear/src/sofa/component/mapping/nonlinear/DistanceMapping.h b/Sofa/Component/Mapping/NonLinear/src/sofa/component/mapping/nonlinear/DistanceMapping.h index aa140874632..595fe8b52b7 100644 --- a/Sofa/Component/Mapping/NonLinear/src/sofa/component/mapping/nonlinear/DistanceMapping.h +++ b/Sofa/Component/Mapping/NonLinear/src/sofa/component/mapping/nonlinear/DistanceMapping.h @@ -23,7 +23,7 @@ #include -#include +#include #include #include #include diff --git a/Sofa/Component/Mapping/NonLinear/src/sofa/component/mapping/nonlinear/DistanceMapping.inl b/Sofa/Component/Mapping/NonLinear/src/sofa/component/mapping/nonlinear/DistanceMapping.inl index 8e2a75fcdfe..a5b7d38ffd6 100644 --- a/Sofa/Component/Mapping/NonLinear/src/sofa/component/mapping/nonlinear/DistanceMapping.inl +++ b/Sofa/Component/Mapping/NonLinear/src/sofa/component/mapping/nonlinear/DistanceMapping.inl @@ -22,7 +22,7 @@ #pragma once #include -#include +#include #include #include #include diff --git a/Sofa/Component/Mapping/NonLinear/src/sofa/component/mapping/nonlinear/SquareDistanceMapping.h b/Sofa/Component/Mapping/NonLinear/src/sofa/component/mapping/nonlinear/SquareDistanceMapping.h index 423e4dea316..138c075f2ea 100644 --- a/Sofa/Component/Mapping/NonLinear/src/sofa/component/mapping/nonlinear/SquareDistanceMapping.h +++ b/Sofa/Component/Mapping/NonLinear/src/sofa/component/mapping/nonlinear/SquareDistanceMapping.h @@ -23,7 +23,7 @@ #include -#include +#include #include #include #include diff --git a/Sofa/Component/Mapping/NonLinear/src/sofa/component/mapping/nonlinear/SquareDistanceMapping.inl b/Sofa/Component/Mapping/NonLinear/src/sofa/component/mapping/nonlinear/SquareDistanceMapping.inl index db6447dd0cb..d016ac706c1 100644 --- a/Sofa/Component/Mapping/NonLinear/src/sofa/component/mapping/nonlinear/SquareDistanceMapping.inl +++ b/Sofa/Component/Mapping/NonLinear/src/sofa/component/mapping/nonlinear/SquareDistanceMapping.inl @@ -22,7 +22,7 @@ #pragma once #include -#include +#include #include #include #include diff --git a/Sofa/Component/Mapping/NonLinear/src/sofa/component/mapping/nonlinear/SquareMapping.h b/Sofa/Component/Mapping/NonLinear/src/sofa/component/mapping/nonlinear/SquareMapping.h index 5c6dc7f701f..a8131b726bc 100644 --- a/Sofa/Component/Mapping/NonLinear/src/sofa/component/mapping/nonlinear/SquareMapping.h +++ b/Sofa/Component/Mapping/NonLinear/src/sofa/component/mapping/nonlinear/SquareMapping.h @@ -23,7 +23,7 @@ #include -#include +#include #include #include diff --git a/Sofa/Component/Mapping/NonLinear/src/sofa/component/mapping/nonlinear/SquareMapping.inl b/Sofa/Component/Mapping/NonLinear/src/sofa/component/mapping/nonlinear/SquareMapping.inl index 46486ce1493..970b25256c2 100644 --- a/Sofa/Component/Mapping/NonLinear/src/sofa/component/mapping/nonlinear/SquareMapping.inl +++ b/Sofa/Component/Mapping/NonLinear/src/sofa/component/mapping/nonlinear/SquareMapping.inl @@ -22,7 +22,7 @@ #pragma once #include -#include +#include #include #include #include diff --git a/Sofa/Component/Mapping/NonLinear/src/sofa/component/mapping/nonlinear/VolumeMapping.h b/Sofa/Component/Mapping/NonLinear/src/sofa/component/mapping/nonlinear/VolumeMapping.h index debd8766eaa..7bb9ba22730 100644 --- a/Sofa/Component/Mapping/NonLinear/src/sofa/component/mapping/nonlinear/VolumeMapping.h +++ b/Sofa/Component/Mapping/NonLinear/src/sofa/component/mapping/nonlinear/VolumeMapping.h @@ -22,7 +22,7 @@ #pragma once #include -#include +#include #include #include #include diff --git a/Sofa/Component/Mapping/NonLinear/src/sofa/component/mapping/nonlinear/VolumeMapping.inl b/Sofa/Component/Mapping/NonLinear/src/sofa/component/mapping/nonlinear/VolumeMapping.inl index c3ce47d3479..899678361e4 100644 --- a/Sofa/Component/Mapping/NonLinear/src/sofa/component/mapping/nonlinear/VolumeMapping.inl +++ b/Sofa/Component/Mapping/NonLinear/src/sofa/component/mapping/nonlinear/VolumeMapping.inl @@ -21,7 +21,7 @@ ******************************************************************************/ #pragma once #include -#include +#include #include namespace sofa::component::mapping::nonlinear From 1b7656adf56c3f847f2de7f599850687e7b8d3fd Mon Sep 17 00:00:00 2001 From: Alex Bilger Date: Thu, 19 Sep 2024 10:17:11 +0200 Subject: [PATCH 12/16] fix --- .../src/sofa/component/mapping/nonlinear/AreaMapping.h | 2 +- .../src/sofa/component/mapping/nonlinear/DistanceMapping.h | 2 +- .../src/sofa/component/mapping/nonlinear/VolumeMapping.h | 2 +- 3 files changed, 3 insertions(+), 3 deletions(-) diff --git a/Sofa/Component/Mapping/NonLinear/src/sofa/component/mapping/nonlinear/AreaMapping.h b/Sofa/Component/Mapping/NonLinear/src/sofa/component/mapping/nonlinear/AreaMapping.h index 2392f310364..01f962962da 100644 --- a/Sofa/Component/Mapping/NonLinear/src/sofa/component/mapping/nonlinear/AreaMapping.h +++ b/Sofa/Component/Mapping/NonLinear/src/sofa/component/mapping/nonlinear/AreaMapping.h @@ -69,7 +69,7 @@ class AreaMapping : public BaseNonLinearMapping const VecCoord_t* m_vertices{nullptr}; - using Inherit1::JacobianEntry; + using JacobianEntry = Inherit1::JacobianEntry; }; #if !defined(SOFA_COMPONENT_MAPPING_NONLINEAR_AREAMAPPING_CPP) diff --git a/Sofa/Component/Mapping/NonLinear/src/sofa/component/mapping/nonlinear/DistanceMapping.h b/Sofa/Component/Mapping/NonLinear/src/sofa/component/mapping/nonlinear/DistanceMapping.h index 595fe8b52b7..f42cfa7eaf2 100644 --- a/Sofa/Component/Mapping/NonLinear/src/sofa/component/mapping/nonlinear/DistanceMapping.h +++ b/Sofa/Component/Mapping/NonLinear/src/sofa/component/mapping/nonlinear/DistanceMapping.h @@ -104,7 +104,7 @@ class DistanceMapping : public BaseNonLinearMapping /// r=b-a only for position (eventual rotation, affine transform... remains null) void computeCoordPositionDifference( Direction& r, const Coord_t& a, const Coord_t& b ); - using Inherit1::JacobianEntry; + using JacobianEntry = Inherit1::JacobianEntry; }; diff --git a/Sofa/Component/Mapping/NonLinear/src/sofa/component/mapping/nonlinear/VolumeMapping.h b/Sofa/Component/Mapping/NonLinear/src/sofa/component/mapping/nonlinear/VolumeMapping.h index 7bb9ba22730..be6e74b0e04 100644 --- a/Sofa/Component/Mapping/NonLinear/src/sofa/component/mapping/nonlinear/VolumeMapping.h +++ b/Sofa/Component/Mapping/NonLinear/src/sofa/component/mapping/nonlinear/VolumeMapping.h @@ -69,7 +69,7 @@ class VolumeMapping : public BaseNonLinearMapping const VecCoord_t* m_vertices{nullptr}; - using Inherit1::JacobianEntry; + using JacobianEntry = Inherit1::JacobianEntry; }; #if !defined(SOFA_COMPONENT_MAPPING_NONLINEAR_VOLUMEMAPPING_CPP) From f67bc84af26ff8591d043a385ec6061c853f8b12 Mon Sep 17 00:00:00 2001 From: Alex Bilger Date: Thu, 19 Sep 2024 10:25:41 +0200 Subject: [PATCH 13/16] second fix --- .../src/sofa/component/mapping/nonlinear/AreaMapping.h | 2 +- .../src/sofa/component/mapping/nonlinear/DistanceMapping.h | 2 +- .../src/sofa/component/mapping/nonlinear/VolumeMapping.h | 2 +- 3 files changed, 3 insertions(+), 3 deletions(-) diff --git a/Sofa/Component/Mapping/NonLinear/src/sofa/component/mapping/nonlinear/AreaMapping.h b/Sofa/Component/Mapping/NonLinear/src/sofa/component/mapping/nonlinear/AreaMapping.h index 01f962962da..3554709407c 100644 --- a/Sofa/Component/Mapping/NonLinear/src/sofa/component/mapping/nonlinear/AreaMapping.h +++ b/Sofa/Component/Mapping/NonLinear/src/sofa/component/mapping/nonlinear/AreaMapping.h @@ -69,7 +69,7 @@ class AreaMapping : public BaseNonLinearMapping const VecCoord_t* m_vertices{nullptr}; - using JacobianEntry = Inherit1::JacobianEntry; + using JacobianEntry = typename Inherit1::JacobianEntry; }; #if !defined(SOFA_COMPONENT_MAPPING_NONLINEAR_AREAMAPPING_CPP) diff --git a/Sofa/Component/Mapping/NonLinear/src/sofa/component/mapping/nonlinear/DistanceMapping.h b/Sofa/Component/Mapping/NonLinear/src/sofa/component/mapping/nonlinear/DistanceMapping.h index f42cfa7eaf2..949fb8fa5d3 100644 --- a/Sofa/Component/Mapping/NonLinear/src/sofa/component/mapping/nonlinear/DistanceMapping.h +++ b/Sofa/Component/Mapping/NonLinear/src/sofa/component/mapping/nonlinear/DistanceMapping.h @@ -104,7 +104,7 @@ class DistanceMapping : public BaseNonLinearMapping /// r=b-a only for position (eventual rotation, affine transform... remains null) void computeCoordPositionDifference( Direction& r, const Coord_t& a, const Coord_t& b ); - using JacobianEntry = Inherit1::JacobianEntry; + using JacobianEntry = typename Inherit1::JacobianEntry; }; diff --git a/Sofa/Component/Mapping/NonLinear/src/sofa/component/mapping/nonlinear/VolumeMapping.h b/Sofa/Component/Mapping/NonLinear/src/sofa/component/mapping/nonlinear/VolumeMapping.h index be6e74b0e04..c8bc2df7b2f 100644 --- a/Sofa/Component/Mapping/NonLinear/src/sofa/component/mapping/nonlinear/VolumeMapping.h +++ b/Sofa/Component/Mapping/NonLinear/src/sofa/component/mapping/nonlinear/VolumeMapping.h @@ -69,7 +69,7 @@ class VolumeMapping : public BaseNonLinearMapping const VecCoord_t* m_vertices{nullptr}; - using JacobianEntry = Inherit1::JacobianEntry; + using JacobianEntry = typename Inherit1::JacobianEntry; }; #if !defined(SOFA_COMPONENT_MAPPING_NONLINEAR_VOLUMEMAPPING_CPP) From 05aa517f6696d70c3413d025dc0273ddb46db9bc Mon Sep 17 00:00:00 2001 From: Alex Bilger Date: Thu, 19 Sep 2024 10:40:19 +0200 Subject: [PATCH 14/16] unused variables --- .../src/sofa/component/mapping/nonlinear/AreaMapping.inl | 2 ++ .../sofa/component/mapping/nonlinear/BaseNonLinearMapping.inl | 4 ++++ .../component/mapping/nonlinear/DistanceFromTargetMapping.inl | 1 + .../component/mapping/nonlinear/SquareDistanceMapping.inl | 1 - .../src/sofa/component/mapping/nonlinear/SquareMapping.inl | 1 - .../src/sofa/component/mapping/nonlinear/VolumeMapping.inl | 2 ++ Sofa/Component/Mapping/NonLinear/tests/AreaMapping_test.cpp | 1 - Sofa/Component/Mapping/NonLinear/tests/VolumeMapping_test.cpp | 1 - 8 files changed, 9 insertions(+), 4 deletions(-) diff --git a/Sofa/Component/Mapping/NonLinear/src/sofa/component/mapping/nonlinear/AreaMapping.inl b/Sofa/Component/Mapping/NonLinear/src/sofa/component/mapping/nonlinear/AreaMapping.inl index 0a343e2fc2e..eed648b526c 100644 --- a/Sofa/Component/Mapping/NonLinear/src/sofa/component/mapping/nonlinear/AreaMapping.inl +++ b/Sofa/Component/Mapping/NonLinear/src/sofa/component/mapping/nonlinear/AreaMapping.inl @@ -137,6 +137,8 @@ template void AreaMapping::apply(const core::MechanicalParams* mparams, DataVecCoord_t& out, const DataVecCoord_t& in) { + SOFA_UNUSED( mparams ); + helper::WriteOnlyAccessor< Data> > _out = out; helper::ReadAccessor< Data > > _in = in; diff --git a/Sofa/Component/Mapping/NonLinear/src/sofa/component/mapping/nonlinear/BaseNonLinearMapping.inl b/Sofa/Component/Mapping/NonLinear/src/sofa/component/mapping/nonlinear/BaseNonLinearMapping.inl index fec4f24a235..b7bdaf79898 100644 --- a/Sofa/Component/Mapping/NonLinear/src/sofa/component/mapping/nonlinear/BaseNonLinearMapping.inl +++ b/Sofa/Component/Mapping/NonLinear/src/sofa/component/mapping/nonlinear/BaseNonLinearMapping.inl @@ -42,6 +42,7 @@ void BaseNonLinearMapping::applyJ( const core::MechanicalParams* mparams, DataVecDeriv_t& out, const DataVecDeriv_t& in) { + SOFA_UNUSED( mparams ); if( jacobian.rowSize() ) { auto dOutWa = sofa::helper::getWriteOnlyAccessor(out); @@ -55,6 +56,7 @@ void BaseNonLinearMapping::applyJT( const core::MechanicalParams* mparams, DataVecDeriv_t& out, const DataVecDeriv_t& in) { + SOFA_UNUSED( mparams ); if( jacobian.rowSize() ) { auto dOutRa = sofa::helper::getReadAccessor(in); @@ -79,6 +81,8 @@ void BaseNonLinearMapping::applyDJT( const core::MechanicalParams* mparams, core::MultiVecDerivId parentForceId, core::ConstMultiVecDerivId childForceId) { + SOFA_UNUSED(childForceId); + const unsigned geometricStiffness = this->d_geometricStiffness.getValue().getSelectedId(); if( !geometricStiffness ) { diff --git a/Sofa/Component/Mapping/NonLinear/src/sofa/component/mapping/nonlinear/DistanceFromTargetMapping.inl b/Sofa/Component/Mapping/NonLinear/src/sofa/component/mapping/nonlinear/DistanceFromTargetMapping.inl index a5fb0255dc6..653f14fb489 100644 --- a/Sofa/Component/Mapping/NonLinear/src/sofa/component/mapping/nonlinear/DistanceFromTargetMapping.inl +++ b/Sofa/Component/Mapping/NonLinear/src/sofa/component/mapping/nonlinear/DistanceFromTargetMapping.inl @@ -184,6 +184,7 @@ void DistanceFromTargetMapping::matrixFreeApplyDJT( const Data>& parentDisplacement, const Data>& childForce) { + SOFA_UNUSED( mparams ); const unsigned geometricStiffness = this->d_geometricStiffness.getValue().getSelectedId(); helper::WriteAccessor parentForceAccessor(parentForce); diff --git a/Sofa/Component/Mapping/NonLinear/src/sofa/component/mapping/nonlinear/SquareDistanceMapping.inl b/Sofa/Component/Mapping/NonLinear/src/sofa/component/mapping/nonlinear/SquareDistanceMapping.inl index d016ac706c1..29361855d0b 100644 --- a/Sofa/Component/Mapping/NonLinear/src/sofa/component/mapping/nonlinear/SquareDistanceMapping.inl +++ b/Sofa/Component/Mapping/NonLinear/src/sofa/component/mapping/nonlinear/SquareDistanceMapping.inl @@ -185,7 +185,6 @@ void SquareDistanceMapping::doUpdateK( const helper::ReadAccessor childForceAccessor(childForce); const SeqEdges& links = l_topology->getEdges(); - unsigned int size = this->fromModel->getSize(); for(size_t i=0; i0) can lead to negative eigen values in geometric stiffness diff --git a/Sofa/Component/Mapping/NonLinear/src/sofa/component/mapping/nonlinear/SquareMapping.inl b/Sofa/Component/Mapping/NonLinear/src/sofa/component/mapping/nonlinear/SquareMapping.inl index 970b25256c2..5433575e8d9 100644 --- a/Sofa/Component/Mapping/NonLinear/src/sofa/component/mapping/nonlinear/SquareMapping.inl +++ b/Sofa/Component/Mapping/NonLinear/src/sofa/component/mapping/nonlinear/SquareMapping.inl @@ -83,7 +83,6 @@ void SquareMapping::doUpdateK(const core::MechanicalParams* mparams, const Data>& childForce, SparseKMatrixEigen& matrix) { SOFA_UNUSED(mparams); - const unsigned geometricStiffness = this->d_geometricStiffness.getValue().getSelectedId(); const helper::ReadAccessor childForceAccessor(childForce); unsigned int size = this->fromModel->getSize(); diff --git a/Sofa/Component/Mapping/NonLinear/src/sofa/component/mapping/nonlinear/VolumeMapping.inl b/Sofa/Component/Mapping/NonLinear/src/sofa/component/mapping/nonlinear/VolumeMapping.inl index 899678361e4..ef2e1bef650 100644 --- a/Sofa/Component/Mapping/NonLinear/src/sofa/component/mapping/nonlinear/VolumeMapping.inl +++ b/Sofa/Component/Mapping/NonLinear/src/sofa/component/mapping/nonlinear/VolumeMapping.inl @@ -116,6 +116,8 @@ template void VolumeMapping::apply(const core::MechanicalParams* mparams, DataVecCoord_t& out, const DataVecCoord_t& in) { + SOFA_UNUSED( mparams ); + helper::WriteOnlyAccessor< Data> > _out = out; helper::ReadAccessor< Data > > _in = in; diff --git a/Sofa/Component/Mapping/NonLinear/tests/AreaMapping_test.cpp b/Sofa/Component/Mapping/NonLinear/tests/AreaMapping_test.cpp index 72e7deb9e5e..1325b49c451 100644 --- a/Sofa/Component/Mapping/NonLinear/tests/AreaMapping_test.cpp +++ b/Sofa/Component/Mapping/NonLinear/tests/AreaMapping_test.cpp @@ -69,7 +69,6 @@ TEST(AreaMapping, firstDerivative) const sofa::type::Mat<3,3,SReal> dA = computeDerivativeArea(vertices); - static constexpr SReal h = 1e-6; for (unsigned int vId = 0; vId < 3; ++vId) { for (unsigned int axis = 0; axis < 3; ++axis) diff --git a/Sofa/Component/Mapping/NonLinear/tests/VolumeMapping_test.cpp b/Sofa/Component/Mapping/NonLinear/tests/VolumeMapping_test.cpp index 26248638920..7a9f7d8b0c7 100644 --- a/Sofa/Component/Mapping/NonLinear/tests/VolumeMapping_test.cpp +++ b/Sofa/Component/Mapping/NonLinear/tests/VolumeMapping_test.cpp @@ -74,7 +74,6 @@ TEST(VolumeMapping, firstDerivative) const sofa::type::Mat<4,3,SReal> dV = computeDerivativeVolume(vertices); - static constexpr SReal h = 1e-6; for (unsigned int vId = 0; vId < 4; ++vId) { for (unsigned int axis = 0; axis < 3; ++axis) From 7ce2e4d633758cc8c7e46a86f590b1df8cd90289 Mon Sep 17 00:00:00 2001 From: Alex Bilger Date: Wed, 25 Sep 2024 13:38:53 +0200 Subject: [PATCH 15/16] rename private members --- .../mapping/nonlinear/BaseNonLinearMapping.h | 4 ++-- .../nonlinear/BaseNonLinearMapping.inl | 20 +++++++++---------- 2 files changed, 12 insertions(+), 12 deletions(-) diff --git a/Sofa/Component/Mapping/NonLinear/src/sofa/component/mapping/nonlinear/BaseNonLinearMapping.h b/Sofa/Component/Mapping/NonLinear/src/sofa/component/mapping/nonlinear/BaseNonLinearMapping.h index 3df90b304ce..822f4595cde 100644 --- a/Sofa/Component/Mapping/NonLinear/src/sofa/component/mapping/nonlinear/BaseNonLinearMapping.h +++ b/Sofa/Component/Mapping/NonLinear/src/sofa/component/mapping/nonlinear/BaseNonLinearMapping.h @@ -90,8 +90,8 @@ class BaseNonLinearMapping : public core::Mapping, public NonLinearMa private: - SparseKMatrixEigen K; ///< Assembled geometric stiffness matrix - type::vector baseMatrices; ///< Jacobian of the mapping, in a vector + SparseKMatrixEigen m_geometricStiffnessMatrix; ///< Assembled geometric stiffness matrix + type::vector m_baseMatrices; }; } diff --git a/Sofa/Component/Mapping/NonLinear/src/sofa/component/mapping/nonlinear/BaseNonLinearMapping.inl b/Sofa/Component/Mapping/NonLinear/src/sofa/component/mapping/nonlinear/BaseNonLinearMapping.inl index b7bdaf79898..93c166d502f 100644 --- a/Sofa/Component/Mapping/NonLinear/src/sofa/component/mapping/nonlinear/BaseNonLinearMapping.inl +++ b/Sofa/Component/Mapping/NonLinear/src/sofa/component/mapping/nonlinear/BaseNonLinearMapping.inl @@ -33,8 +33,8 @@ init() { core::Mapping::init(); - this->baseMatrices.resize( 1 ); - this->baseMatrices[0] = &this->jacobian; + m_baseMatrices.resize( 1 ); + m_baseMatrices[0] = &this->jacobian; } template @@ -94,11 +94,11 @@ void BaseNonLinearMapping::applyDJT( const SReal kFactor = mparams->kFactor(); - if( K.compressedMatrix.nonZeros() ) + if( m_geometricStiffnessMatrix.compressedMatrix.nonZeros() ) { helper::WriteAccessor parentForceAccessor(parentForce); helper::ReadAccessor parentDisplacementAccessor(parentDisplacement); - K.addMult( parentForceAccessor.wref(), parentDisplacementAccessor.ref(), static_cast(kFactor) ); + m_geometricStiffnessMatrix.addMult( parentForceAccessor.wref(), parentDisplacementAccessor.ref(), static_cast(kFactor) ); } else { @@ -112,13 +112,13 @@ void BaseNonLinearMapping::applyDJT( template const linearalgebra::BaseMatrix* BaseNonLinearMapping::getK() { - return &K; + return &m_geometricStiffnessMatrix; } template const type::vector* BaseNonLinearMapping::getJs() { - return &baseMatrices; + return &m_baseMatrices; } template @@ -126,18 +126,18 @@ void BaseNonLinearMapping:: updateK(const core::MechanicalParams* mparams, core::ConstMultiVecDerivId childForceId) { const unsigned geometricStiffness = this->d_geometricStiffness.getValue().getSelectedId(); - if( !geometricStiffness ) { this->K.resize(0,0); return; } + if( !geometricStiffness ) { this->m_geometricStiffnessMatrix.resize(0,0); return; } const Data >& childForce = *childForceId[this->toModel.get()].read(); { unsigned int kSize = this->fromModel->getSize(); - K.resizeBlocks(kSize, kSize); + m_geometricStiffnessMatrix.resizeBlocks(kSize, kSize); } - doUpdateK(mparams, childForce, K); + doUpdateK(mparams, childForce, m_geometricStiffnessMatrix); - K.compress(); + m_geometricStiffnessMatrix.compress(); } From 37649f95e53fb0a69e77b2f0fa8b0a0481bbe50d Mon Sep 17 00:00:00 2001 From: Alex Bilger Date: Mon, 30 Sep 2024 14:13:22 +0200 Subject: [PATCH 16/16] change protected member --- .../component/mapping/nonlinear/AreaMapping.inl | 10 +++++----- .../mapping/nonlinear/BaseNonLinearMapping.h | 2 +- .../mapping/nonlinear/BaseNonLinearMapping.inl | 12 ++++++------ .../nonlinear/DistanceFromTargetMapping.inl | 8 ++++---- .../mapping/nonlinear/DistanceMapping.inl | 10 +++++----- .../mapping/nonlinear/SquareDistanceMapping.inl | 14 +++++++------- .../component/mapping/nonlinear/SquareMapping.inl | 10 +++++----- .../component/mapping/nonlinear/VolumeMapping.inl | 10 +++++----- 8 files changed, 38 insertions(+), 38 deletions(-) diff --git a/Sofa/Component/Mapping/NonLinear/src/sofa/component/mapping/nonlinear/AreaMapping.inl b/Sofa/Component/Mapping/NonLinear/src/sofa/component/mapping/nonlinear/AreaMapping.inl index eed648b526c..5b856d8b7ee 100644 --- a/Sofa/Component/Mapping/NonLinear/src/sofa/component/mapping/nonlinear/AreaMapping.inl +++ b/Sofa/Component/Mapping/NonLinear/src/sofa/component/mapping/nonlinear/AreaMapping.inl @@ -123,7 +123,7 @@ void AreaMapping::init() typename core::behavior::MechanicalState::ReadVecCoord pos = this->getFromModel()->readPositions(); this->getToModel()->resize( nbTriangles ); - this->jacobian.resizeBlocks(nbTriangles, pos.size()); + this->m_jacobian.resizeBlocks(nbTriangles, pos.size()); Inherit1::init(); @@ -146,7 +146,7 @@ void AreaMapping::apply(const core::MechanicalParams* mparams, const auto& triangles = l_topology->getTriangles(); - this->jacobian.clear(); + this->m_jacobian.clear(); for (unsigned int triangleId = 0; triangleId < triangles.size(); ++triangleId) { @@ -176,17 +176,17 @@ void AreaMapping::apply(const core::MechanicalParams* mparams, //insertion in increasing column order std::sort(jacobianEntries.begin(), jacobianEntries.end()); - this->jacobian.beginRow(triangleId); + this->m_jacobian.beginRow(triangleId); for (const auto& [vertexId, jacobianValue] : jacobianEntries) { for (unsigned d = 0; d < In::spatial_dimensions; ++d) { - this->jacobian.insertBack(triangleId, vertexId * Nin + d, jacobianValue[d]); + this->m_jacobian.insertBack(triangleId, vertexId * Nin + d, jacobianValue[d]); } } } - this->jacobian.compress(); + this->m_jacobian.compress(); } template diff --git a/Sofa/Component/Mapping/NonLinear/src/sofa/component/mapping/nonlinear/BaseNonLinearMapping.h b/Sofa/Component/Mapping/NonLinear/src/sofa/component/mapping/nonlinear/BaseNonLinearMapping.h index 822f4595cde..847c8628cdf 100644 --- a/Sofa/Component/Mapping/NonLinear/src/sofa/component/mapping/nonlinear/BaseNonLinearMapping.h +++ b/Sofa/Component/Mapping/NonLinear/src/sofa/component/mapping/nonlinear/BaseNonLinearMapping.h @@ -67,7 +67,7 @@ class BaseNonLinearMapping : public core::Mapping, public NonLinearMa const Data >& parentDisplacement, const Data >& childForce) = 0; - SparseMatrixEigen jacobian; ///< Jacobian of the mapping + SparseMatrixEigen m_jacobian; ///< Jacobian of the mapping virtual void doUpdateK( const core::MechanicalParams* mparams, const Data >& childForce, diff --git a/Sofa/Component/Mapping/NonLinear/src/sofa/component/mapping/nonlinear/BaseNonLinearMapping.inl b/Sofa/Component/Mapping/NonLinear/src/sofa/component/mapping/nonlinear/BaseNonLinearMapping.inl index 93c166d502f..742d3303e87 100644 --- a/Sofa/Component/Mapping/NonLinear/src/sofa/component/mapping/nonlinear/BaseNonLinearMapping.inl +++ b/Sofa/Component/Mapping/NonLinear/src/sofa/component/mapping/nonlinear/BaseNonLinearMapping.inl @@ -34,7 +34,7 @@ init() core::Mapping::init(); m_baseMatrices.resize( 1 ); - m_baseMatrices[0] = &this->jacobian; + m_baseMatrices[0] = &this->m_jacobian; } template @@ -43,11 +43,11 @@ void BaseNonLinearMapping::applyJ( const DataVecDeriv_t& in) { SOFA_UNUSED( mparams ); - if( jacobian.rowSize() ) + if( m_jacobian.rowSize() ) { auto dOutWa = sofa::helper::getWriteOnlyAccessor(out); auto dInRa = sofa::helper::getReadAccessor(in); - jacobian.mult(dOutWa.wref(),dInRa.ref()); + m_jacobian.mult(dOutWa.wref(),dInRa.ref()); } } @@ -57,11 +57,11 @@ void BaseNonLinearMapping::applyJT( const DataVecDeriv_t& in) { SOFA_UNUSED( mparams ); - if( jacobian.rowSize() ) + if( m_jacobian.rowSize() ) { auto dOutRa = sofa::helper::getReadAccessor(in); auto dInWa = sofa::helper::getWriteOnlyAccessor(out); - jacobian.addMultTranspose(dInWa.wref(),dOutRa.ref()); + m_jacobian.addMultTranspose(dInWa.wref(),dOutRa.ref()); } } @@ -73,7 +73,7 @@ void BaseNonLinearMapping::applyJT( SOFA_UNUSED(cparams); auto childMatRa = sofa::helper::getReadAccessor(in); auto parentMatWa = sofa::helper::getWriteAccessor(out); - addMultTransposeEigen(parentMatWa.wref(), jacobian.compressedMatrix, childMatRa.ref()); + addMultTransposeEigen(parentMatWa.wref(), m_jacobian.compressedMatrix, childMatRa.ref()); } template diff --git a/Sofa/Component/Mapping/NonLinear/src/sofa/component/mapping/nonlinear/DistanceFromTargetMapping.inl b/Sofa/Component/Mapping/NonLinear/src/sofa/component/mapping/nonlinear/DistanceFromTargetMapping.inl index 653f14fb489..c350fdf6a0d 100644 --- a/Sofa/Component/Mapping/NonLinear/src/sofa/component/mapping/nonlinear/DistanceFromTargetMapping.inl +++ b/Sofa/Component/Mapping/NonLinear/src/sofa/component/mapping/nonlinear/DistanceFromTargetMapping.inl @@ -137,7 +137,7 @@ void DistanceFromTargetMapping::apply(const core::MechanicalParams * const helper::ReadAccessor indices(d_indices); const helper::ReadAccessor targetPositions(d_targetPositions); - this->jacobian.resizeBlocks(out.size(),in.size()); + this->m_jacobian.resizeBlocks(out.size(),in.size()); directions.resize(out.size()); invlengths.resize(out.size()); @@ -166,15 +166,15 @@ void DistanceFromTargetMapping::apply(const core::MechanicalParams * for(unsigned j=0; jjacobian.beginRow(i*Nout+j); + this->m_jacobian.beginRow(i*Nout+j); for(unsigned k=0; kjacobian.insertBack( i*Nout+j, indices[i]*Nin+k, gap[k] ); + this->m_jacobian.insertBack( i*Nout+j, indices[i]*Nin+k, gap[k] ); } } } - this->jacobian.compress(); + this->m_jacobian.compress(); } template diff --git a/Sofa/Component/Mapping/NonLinear/src/sofa/component/mapping/nonlinear/DistanceMapping.inl b/Sofa/Component/Mapping/NonLinear/src/sofa/component/mapping/nonlinear/DistanceMapping.inl index a5b7d38ffd6..6c90d4616fe 100644 --- a/Sofa/Component/Mapping/NonLinear/src/sofa/component/mapping/nonlinear/DistanceMapping.inl +++ b/Sofa/Component/Mapping/NonLinear/src/sofa/component/mapping/nonlinear/DistanceMapping.inl @@ -70,7 +70,7 @@ void DistanceMapping::init() typename core::behavior::MechanicalState::ReadVecCoord pos = this->getFromModel()->readPositions(); this->getToModel()->resize( links.size() ); - this->jacobian.resizeBlocks(links.size(),pos.size()); + this->m_jacobian.resizeBlocks(links.size(),pos.size()); directions.resize(links.size()); invlengths.resize(links.size()); @@ -108,7 +108,7 @@ void DistanceMapping::apply(const core::MechanicalParams * /*mparams* helper::ReadAccessor > > restLengths(d_restLengths); const SeqEdges& links = l_topology->getEdges(); - this->jacobian.clear(); + this->m_jacobian.clear(); for (unsigned int i = 0; i < links.size(); ++i) { @@ -144,17 +144,17 @@ void DistanceMapping::apply(const core::MechanicalParams * /*mparams* //invert to insert in increasing column order std::sort(jacobianEntries.begin(), jacobianEntries.end()); - this->jacobian.beginRow(i); + this->m_jacobian.beginRow(i); for (const auto& [vertexId, jacobianValue] : jacobianEntries) { for (unsigned k = 0; k < In::spatial_dimensions; ++k) { - this->jacobian.insertBack(i, vertexId * Nin + k, jacobianValue[k]); + this->m_jacobian.insertBack(i, vertexId * Nin + k, jacobianValue[k]); } } } - this->jacobian.compress(); + this->m_jacobian.compress(); } template diff --git a/Sofa/Component/Mapping/NonLinear/src/sofa/component/mapping/nonlinear/SquareDistanceMapping.inl b/Sofa/Component/Mapping/NonLinear/src/sofa/component/mapping/nonlinear/SquareDistanceMapping.inl index 29361855d0b..54d11d5cb55 100644 --- a/Sofa/Component/Mapping/NonLinear/src/sofa/component/mapping/nonlinear/SquareDistanceMapping.inl +++ b/Sofa/Component/Mapping/NonLinear/src/sofa/component/mapping/nonlinear/SquareDistanceMapping.inl @@ -88,7 +88,7 @@ void SquareDistanceMapping::apply( const helper::ReadAccessor in (dIn); const SeqEdges& links = l_topology->getEdges(); - this->jacobian.resizeBlocks(out.size(),in.size()); + this->m_jacobian.resizeBlocks(out.size(),in.size()); Direction gap; @@ -118,24 +118,24 @@ void SquareDistanceMapping::apply( // } - this->jacobian.beginRow(i); + this->m_jacobian.beginRow(i); if( links[i][1]jacobian.insertBack( i, links[i][1]*Nin+k, gap[k] ); + this->m_jacobian.insertBack( i, links[i][1]*Nin+k, gap[k] ); for(unsigned k=0; kjacobian.insertBack( i, links[i][0]*Nin+k, -gap[k] ); + this->m_jacobian.insertBack( i, links[i][0]*Nin+k, -gap[k] ); } else { for(unsigned k=0; kjacobian.insertBack( i, links[i][0]*Nin+k, -gap[k] ); + this->m_jacobian.insertBack( i, links[i][0]*Nin+k, -gap[k] ); for(unsigned k=0; kjacobian.insertBack( i, links[i][1]*Nin+k, gap[k] ); + this->m_jacobian.insertBack( i, links[i][1]*Nin+k, gap[k] ); } } - this->jacobian.compress(); + this->m_jacobian.compress(); } template diff --git a/Sofa/Component/Mapping/NonLinear/src/sofa/component/mapping/nonlinear/SquareMapping.inl b/Sofa/Component/Mapping/NonLinear/src/sofa/component/mapping/nonlinear/SquareMapping.inl index 5433575e8d9..a74203abe0c 100644 --- a/Sofa/Component/Mapping/NonLinear/src/sofa/component/mapping/nonlinear/SquareMapping.inl +++ b/Sofa/Component/Mapping/NonLinear/src/sofa/component/mapping/nonlinear/SquareMapping.inl @@ -42,19 +42,19 @@ void SquareMapping::apply(const core::MechanicalParams* mparams, size_t size = in.size(); this->getToModel()->resize( size ); - this->jacobian.resizeBlocks( size, size ); - this->jacobian.reserve( size ); + this->m_jacobian.resizeBlocks( size, size ); + this->m_jacobian.reserve( size ); for( unsigned i=0 ; ijacobian.beginRow(i); - this->jacobian.insertBack( i, i, 2.0*x ); + this->m_jacobian.beginRow(i); + this->m_jacobian.insertBack( i, i, 2.0*x ); } - this->jacobian.compress(); + this->m_jacobian.compress(); } template diff --git a/Sofa/Component/Mapping/NonLinear/src/sofa/component/mapping/nonlinear/VolumeMapping.inl b/Sofa/Component/Mapping/NonLinear/src/sofa/component/mapping/nonlinear/VolumeMapping.inl index ef2e1bef650..06d97729f97 100644 --- a/Sofa/Component/Mapping/NonLinear/src/sofa/component/mapping/nonlinear/VolumeMapping.inl +++ b/Sofa/Component/Mapping/NonLinear/src/sofa/component/mapping/nonlinear/VolumeMapping.inl @@ -102,7 +102,7 @@ void VolumeMapping::init() typename core::behavior::MechanicalState::ReadVecCoord pos = this->getFromModel()->readPositions(); this->getToModel()->resize( nbTetrahedra ); - this->jacobian.resizeBlocks(nbTetrahedra, pos.size()); + this->m_jacobian.resizeBlocks(nbTetrahedra, pos.size()); Inherit1::init(); @@ -125,7 +125,7 @@ void VolumeMapping::apply(const core::MechanicalParams* mparams, const auto& tetrahedra = l_topology->getTetrahedra(); - this->jacobian.clear(); + this->m_jacobian.clear(); for (unsigned int tetId = 0; tetId < tetrahedra.size(); ++tetId) { @@ -165,17 +165,17 @@ void VolumeMapping::apply(const core::MechanicalParams* mparams, //insertion in increasing column order std::sort(jacobianEntries.begin(), jacobianEntries.end()); - this->jacobian.beginRow(tetId); + this->m_jacobian.beginRow(tetId); for (const auto& [vertexId, jacobianValue] : jacobianEntries) { for (unsigned d = 0; d < In::spatial_dimensions; ++d) { - this->jacobian.insertBack(tetId, vertexId * Nin + d, jacobianValue[d]); + this->m_jacobian.insertBack(tetId, vertexId * Nin + d, jacobianValue[d]); } } } - this->jacobian.compress(); + this->m_jacobian.compress(); } template