Sofa vDot product implementations needs overloading operator*

[jjcasmar]

When adding a new DataType for a MechanicalObject, I need to provide an operator* implementations which computes the dot product between two DataTypes.

    Real r = 0.0;

    if (a.type == sofa::core::V_COORD && b.type == sofa::core::V_COORD)
    {
        const VecCoord &va = this->read(core::ConstVecCoordId(a))->getValue();
        const VecCoord &vb = this->read(core::ConstVecCoordId(b))->getValue();

        for (unsigned int i=0; i<va.size(); i++)
        {
            r += va[i] * vb[i];
        }
    }

This is clearly an issue if for example you use Eigen as a DataType, as the operator has a clear meaning. If for example you have an Eigen::Matrix<double, 3, 3>, you cant override operator because you wont be able to perform normal algebraic matrix multiplication.

Consider using free functions instead of overloading operator*.

Real r = 0.0;

if (a.type == sofa::core::V_COORD && b.type == sofa::core::V_COORD)
{
    const VecCoord &va = this->read(core::ConstVecCoordId(a))->getValue();
    const VecCoord &vb = this->read(core::ConstVecCoordId(b))->getValue();
    using Sofa::dot;
    for (unsigned int i=0; i<va.size(); i++)
    {
        r += dot(va[i], vb[i])
    }
}

[jnbrunet]

Hey @jjcasmar

Thank you, this is indeed a good suggestion. There are probably other places where such operator overload is used and should be replaced by a free function. We should decide on where to place these free functions and maybe try to list all operator overloads usage in the rest of SOFA and its plugins.

J-N

[jjcasmar]

Probably the best way to see where this operator overloads are used is to remove them from StdVectorTypes (Vec3Types and family) and see where we get the compile errors.

[jjcasmar]

For the record, what I am doing to alleviate this issue is to specialize the MechanicalObject functions that fail to compile with the default implementation and my DataType. There are a couple that fail because the access method operator[] doesn't make sense or can't be overriden (operator[] exists for Eigen::Matrix so its ambiguos for example). This is the list I had to specialize

void MechanicalObject<VNCS::F33D>::vThreshold(sofa::core::VecId v, double t)
{
    // Do nothing
}
template <>
double MechanicalObject<VNCS::F33D>::vDot(const core::ExecParams *, core::ConstVecId a, core::ConstVecId b)
{
    Real r = 0.0;
    return r;
}

template <>
double MechanicalObject<VNCS::F33D>::getConstraintJacobianTimesVecDeriv(unsigned int line, core::ConstVecId id)
{
    SReal result = 0;
    return result;
}

template <>
void MechanicalObject<VNCS::F33D>::applyScale(const double sx, const double sy, const double sz)
{
}
template <>
void MechanicalObject<VNCS::F33D>::storeResetState()
{
    // store a reset state only for independent dofs (mapped dofs are deduced from independent dofs)
    if (!isIndependent())
        return;
}

template <>
void MechanicalObject<VNCS::F33D>::buildIdentityBlocksInJacobian(const sofa::helper::vector<unsigned int> &list_n,
                                                                 core::MatrixDerivId &mID)
{
}

template <>
SReal MechanicalObject<VNCS::F33D>::vSum(const core::ExecParams *params, core::ConstVecId a, unsigned l)
{
    Real r = 0.0;
    return r;
}

template <>
SReal MechanicalObject<VNCS::F33D>::vMax(const core::ExecParams *params, core::ConstVecId a)
{
    Real r = 0.0;
    return r;
}

template <>
void MechanicalObject<VNCS::F33D>::getConstraintJacobian(const core::ConstraintParams *cParams,
                                                         sofa::defaulttype::BaseMatrix *J,
                                                         unsigned int &off)
{
}

Im not sure if this specializations make sense, but I think they do as in my particular case VNCS::F33D is a dependent DoF and I am using a CG solver, so I think most of them are not called anyway.

I found that if a DataType can only be used as a dependent DoF (for example, in this case, VNCS::F33D is a deformation gradient which is by definition a dependent DoF), that storeResetState could be specialized.