CFD/eigen
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

Cross product for vectors of size 2. Fixes #1037

Browse Source
This commit is contained in:
Gabriele Buondonno 2022-11-15 22:39:42 +00:00 committed by Antonio Sánchez
parent 8588d8c74b
commit 6431dfdb50
7 changed files with 108 additions and 34 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

13
Eigen/src/Core/MatrixBase.h
View File

@ -387,20 +387,9 @@ template<typename Derived> class MatrixBase
/////////// Geometry module ///////////
#ifndef EIGEN_PARSED_BY_DOXYGEN
/// \internal helper struct to form the return type of the cross product
template<typename OtherDerived> struct cross_product_return_type {
typedef typename ScalarBinaryOpTraits<typename internal::traits<Derived>::Scalar,typename internal::traits<OtherDerived>::Scalar>::ReturnType Scalar;
typedef Matrix<Scalar,MatrixBase::RowsAtCompileTime,MatrixBase::ColsAtCompileTime> type;
};
#endif // EIGEN_PARSED_BY_DOXYGEN
template<typename OtherDerived>
EIGEN_DEVICE_FUNC
#ifndef EIGEN_PARSED_BY_DOXYGEN
inline typename cross_product_return_type<OtherDerived>::type
#else
inline PlainObject
#endif
inline typename internal::cross_impl<Derived, OtherDerived>::return_type
cross(const MatrixBase<OtherDerived>& other) const;
template<typename OtherDerived>

4
Eigen/src/Core/util/ForwardDeclarations.h
View File

@ -270,8 +270,10 @@ template<typename VectorsType, typename CoeffsType, int Side=OnTheLeft> class Ho
template<typename Scalar> class JacobiRotation;
// Geometry module:
namespace internal {
template<typename Derived, typename OtherDerived, int Size = MatrixBase<Derived>::SizeAtCompileTime> struct cross_impl;
}
template<typename Derived, int Dim_> class RotationBase;
template<typename Lhs, typename Rhs> class Cross;
template<typename Derived> class QuaternionBase;
template<typename Scalar> class Rotation2D;
template<typename Scalar> class AngleAxis;

82
Eigen/src/Geometry/OrthoMethods.h
View File

@ -15,39 +15,83 @@
namespace Eigen {
namespace internal {
// Vector3 version (default)
template<typename Derived, typename OtherDerived, int Size>
struct cross_impl
{
typedef typename ScalarBinaryOpTraits<typename internal::traits<Derived>::Scalar,typename internal::traits<OtherDerived>::Scalar>::ReturnType Scalar;
typedef Matrix<Scalar,MatrixBase<Derived>::RowsAtCompileTime,MatrixBase<Derived>::ColsAtCompileTime> return_type;
static EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE
return_type run(const MatrixBase<Derived>& first, const MatrixBase<OtherDerived>& second)
{
EIGEN_STATIC_ASSERT_VECTOR_SPECIFIC_SIZE(Derived,3)
EIGEN_STATIC_ASSERT_VECTOR_SPECIFIC_SIZE(OtherDerived,3)
// Note that there is no need for an expression here since the compiler
// optimize such a small temporary very well (even within a complex expression)
typename internal::nested_eval<Derived,2>::type lhs(first.derived());
typename internal::nested_eval<OtherDerived,2>::type rhs(second.derived());
return return_type(
numext::conj(lhs.coeff(1) * rhs.coeff(2) - lhs.coeff(2) * rhs.coeff(1)),
numext::conj(lhs.coeff(2) * rhs.coeff(0) - lhs.coeff(0) * rhs.coeff(2)),
numext::conj(lhs.coeff(0) * rhs.coeff(1) - lhs.coeff(1) * rhs.coeff(0))
);
}
};
// Vector2 version
template<typename Derived, typename OtherDerived>
struct cross_impl<Derived, OtherDerived, 2>
{
typedef typename ScalarBinaryOpTraits<typename internal::traits<Derived>::Scalar,typename internal::traits<OtherDerived>::Scalar>::ReturnType Scalar;
typedef Scalar return_type;
static EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE
return_type run(const MatrixBase<Derived>& first, const MatrixBase<OtherDerived>& second)
{
EIGEN_STATIC_ASSERT_VECTOR_SPECIFIC_SIZE(Derived,2);
EIGEN_STATIC_ASSERT_VECTOR_SPECIFIC_SIZE(OtherDerived,2);
typename internal::nested_eval<Derived,2>::type lhs(first.derived());
typename internal::nested_eval<OtherDerived,2>::type rhs(second.derived());
return numext::conj(lhs.coeff(0) * rhs.coeff(1) - lhs.coeff(1) * rhs.coeff(0));
}
};
} // end namespace internal
/** \geometry_module \ingroup Geometry_Module
*
* \returns the cross product of \c *this and \a other
* \returns the cross product of \c *this and \a other. This is either a scalar for size-2 vectors or a size-3 vector for size-3 vectors.
*
* Here is a very good explanation of cross-product: http://xkcd.com/199/
* This method is implemented for two different cases: between vectors of fixed size 2 and between vectors of fixed size 3.
*
* With complex numbers, the cross product is implemented as
* \f$ (\mathbf{a}+i\mathbf{b}) \times (\mathbf{c}+i\mathbf{d}) = (\mathbf{a} \times \mathbf{c} - \mathbf{b} \times \mathbf{d}) - i(\mathbf{a} \times \mathbf{d} - \mathbf{b} \times \mathbf{c})\f$
* For vectors of size 3, the output is simply the traditional cross product.
*
* For vectors of size 2, the output is a scalar.
* Given vectors \f$ v = \begin{bmatrix} v_1 & v_2 \end{bmatrix} \f$ and \f$ w = \begin{bmatrix} w_1 & w_2 \end{bmatrix} \f$,
* the result is simply \f$ v\times w = \overline{v_1 w_2 - v_2 w_1} = \text{conj}\left|\begin{smallmatrix} v_1 & w_1 \\ v_2 & w_2 \end{smallmatrix}\right| \f$;
* or, to put it differently, it is the third coordinate of the cross product of \f$ \begin{bmatrix} v_1 & v_2 & v_3 \end{bmatrix} \f$ and \f$ \begin{bmatrix} w_1 & w_2 & w_3 \end{bmatrix} \f$.
* For real-valued inputs, the result can be interpreted as the signed area of a parallelogram spanned by the two vectors.
*
* \note With complex numbers, the cross product is implemented as
* \f$ (\mathbf{a}+i\mathbf{b}) \times (\mathbf{c}+i\mathbf{d}) = (\mathbf{a} \times \mathbf{c} - \mathbf{b} \times \mathbf{d}) - i(\mathbf{a} \times \mathbf{d} + \mathbf{b} \times \mathbf{c})\f$
*
* \sa MatrixBase::cross3()
*/
template<typename Derived>
template<typename OtherDerived>
#ifndef EIGEN_PARSED_BY_DOXYGEN
EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE
typename MatrixBase<Derived>::template cross_product_return_type<OtherDerived>::type
#ifndef EIGEN_PARSED_BY_DOXYGEN
typename internal::cross_impl<Derived, OtherDerived>::return_type
#else
typename MatrixBase<Derived>::PlainObject
inline std::conditional_t<SizeAtCompileTime==2, Scalar, PlainObject>
#endif
MatrixBase<Derived>::cross(const MatrixBase<OtherDerived>& other) const
{
EIGEN_STATIC_ASSERT_VECTOR_SPECIFIC_SIZE(Derived,3)
EIGEN_STATIC_ASSERT_VECTOR_SPECIFIC_SIZE(OtherDerived,3)
// Note that there is no need for an expression here since the compiler
// optimize such a small temporary very well (even within a complex expression)
typename internal::nested_eval<Derived,2>::type lhs(derived());
typename internal::nested_eval<OtherDerived,2>::type rhs(other.derived());
return typename cross_product_return_type<OtherDerived>::type(
numext::conj(lhs.coeff(1) * rhs.coeff(2) - lhs.coeff(2) * rhs.coeff(1)),
numext::conj(lhs.coeff(2) * rhs.coeff(0) - lhs.coeff(0) * rhs.coeff(2)),
numext::conj(lhs.coeff(0) * rhs.coeff(1) - lhs.coeff(1) * rhs.coeff(0))
);
return internal::cross_impl<Derived, OtherDerived>::run(*this, other);
}
namespace internal {

3
doc/QuickReference.dox
View File

@ -367,7 +367,8 @@ vec2 = vec1.normalized(); vec1.normalize(); // inplace \endcode
<tr class="alt"><td>
\link MatrixBase::cross() cross product \endlink \matrixworld</td><td>\code
#include <Eigen/Geometry>
vec3 = vec1.cross(vec2);\endcode</td></tr>
v3c = v3a.cross(v3b); // size-3 vectors
scalar = v2a.cross(v2b); // size-2 vectors \endcode</td></tr>
</table>
<a href="#" class="top">top</a>

2
doc/TutorialMatrixArithmetic.dox
View File

@ -158,7 +158,7 @@ For dot product and cross product, you need the \link MatrixBase::dot() dot()\en
\verbinclude tut_arithmetic_dot_cross.out
</td></tr></table>
Remember that cross product is only for vectors of size 3. Dot product is for vectors of any sizes.
Cross product is defined in Eigen not only for vectors of size 3 but also for those of size 2, check \link MatrixBase::cross() the doc\endlink for details. Dot product is for vectors of any sizes.
When using complex numbers, Eigen's dot product is conjugate-linear in the first variable and linear in the
second variable.

5
doc/examples/tut_arithmetic_dot_cross.cpp
View File

@ -9,5 +9,10 @@ int main()
std::cout << "Dot product: " << v.dot(w) << std::endl;
double dp = v.adjoint()*w; // automatic conversion of the inner product to a scalar
std::cout << "Dot product via a matrix product: " << dp << std::endl;
std::cout << "Cross product:\n" << v.cross(w) << std::endl;
Eigen::Vector2d v2(1,2);
Eigen::Vector2d w2(0,1);
double cp = v2.cross(w2); // returning a scalar between size-2 vectors
std::cout << "Cross product for 2D vectors: " << cp << std::endl;
}

33
test/geo_orthomethods.cpp
View File

@ -78,6 +78,36 @@ template<typename Scalar> void orthomethods_3()
VERIFY_IS_APPROX(v2.cross(v1), rv1.cross(v1));
}
template<typename Scalar> void orthomethods_2()
{
typedef typename NumTraits<Scalar>::Real RealScalar;
typedef Matrix<Scalar,2,1> Vector2;
typedef Matrix<Scalar,3,1> Vector3;
Vector3 v30 = Vector3::Random(),
v31 = Vector3::Random();
Vector2 v20 = v30.template head<2>();
Vector2 v21 = v31.template head<2>();
VERIFY_IS_MUCH_SMALLER_THAN(v20.cross(v20), Scalar(1));
VERIFY_IS_MUCH_SMALLER_THAN(v21.cross(v21), Scalar(1));
VERIFY_IS_APPROX(v20.cross(v21), v30.cross(v31).z());
Vector2 v20Rot90(numext::conj(-v20.y()), numext::conj(v20.x()));
VERIFY_IS_APPROX(v20.cross( v20Rot90), v20.squaredNorm());
VERIFY_IS_APPROX(v20.cross(-v20Rot90), -v20.squaredNorm());
Vector2 v21Rot90(numext::conj(-v21.y()), numext::conj(v21.x()));
VERIFY_IS_APPROX(v21.cross( v21Rot90), v21.squaredNorm());
VERIFY_IS_APPROX(v21.cross(-v21Rot90), -v21.squaredNorm());
// check mixed product
typedef Matrix<RealScalar, 2, 1> RealVector2;
RealVector2 rv21 = RealVector2::Random();
v21 = rv21.template cast<Scalar>();
VERIFY_IS_APPROX(v20.cross(v21), v20.cross(rv21));
VERIFY_IS_APPROX(v21.cross(v20), rv21.cross(v20));
}
template<typename Scalar, int Size> void orthomethods(int size=Size)
{
typedef typename NumTraits<Scalar>::Real RealScalar;
@ -119,6 +149,9 @@ template<typename Scalar, int Size> void orthomethods(int size=Size)
EIGEN_DECLARE_TEST(geo_orthomethods)
{
for(int i = 0; i < g_repeat; i++) {
CALL_SUBTEST_1( orthomethods_2<float>() );
CALL_SUBTEST_2( orthomethods_2<double>() );
CALL_SUBTEST_4( orthomethods_2<std::complex<double> >() );
CALL_SUBTEST_1( orthomethods_3<float>() );
CALL_SUBTEST_2( orthomethods_3<double>() );
CALL_SUBTEST_4( orthomethods_3<std::complex<double> >() );