From ede066f8522657ac141c46b9e5a6c920ccc2a52b Mon Sep 17 00:00:00 2001 From: PabloIbannez Date: Mon, 25 Dec 2023 10:40:40 +0100 Subject: [PATCH] Added qualifiers --- src/utils/quaternion.cuh | 196 ++++++++++++++++++++------------------- 1 file changed, 100 insertions(+), 96 deletions(-) diff --git a/src/utils/quaternion.cuh b/src/utils/quaternion.cuh index 2c384549..bdaec2a4 100644 --- a/src/utils/quaternion.cuh +++ b/src/utils/quaternion.cuh @@ -1,7 +1,7 @@ /* P. Palacios Alonso 2021 Some quaternion algebra and useful functions Notation: - n(real) - Scalar part + n(real) - Scalar part v(real3) - vectorial part q(real4) - quaternion */ @@ -13,146 +13,150 @@ Notation: namespace uammd{ class Quat{ - public: - - real n; - real3 v; - - QUATTR Quat(); - QUATTR Quat(real4 q); - QUATTR Quat(real n, real3 v); - QUATTR Quat(real n, real vx, real vy, real vz); - - QUATTR Quat operator+(Quat q1); - QUATTR void operator+=(Quat q); - QUATTR Quat operator-(Quat q); - QUATTR void operator-=(Quat q); - QUATTR Quat operator*(Quat q); - QUATTR void operator*=(Quat q); - QUATTR Quat operator*(real scalar); - QUATTR void operator*=(real scalar); - QUATTR Quat operator/(real scalar); - QUATTR void operator/=(real scalar); - QUATTR void operator=(real4 q); - - QUATTR real3 getVx(); //Returns the first vector of the reference system encoded by the quaternion - QUATTR real3 getVy(); //Returns the second vector of the reference system encoded by the quaternion - QUATTR real3 getVz(); //Returns the third vector of the reference system encoded by the quaternion - QUATTR real4 to_real4(); - - QUATTR Quat getConjugate(); + public: + + real n; + real3 v; + + QUATTR Quat(); + QUATTR Quat(const real4& q); + QUATTR Quat(real n, const real3& v); + QUATTR Quat(real n, real vx, real vy, real vz); + + QUATTR Quat operator+ (const Quat& q) const; + QUATTR void operator+= (const Quat& q); + QUATTR Quat operator- (const Quat& q) const; + QUATTR void operator-= (const Quat& q); + QUATTR Quat operator* (const Quat& q) const; + QUATTR void operator*= (const Quat& q); + QUATTR Quat operator* (real scalar) const; + QUATTR void operator*= (real scalar); + QUATTR Quat operator/ (real scalar) const; + QUATTR void operator/= (real scalar); + QUATTR void operator= (const real4& q); + + QUATTR real3 getVx() const; + QUATTR real3 getVy() const; + QUATTR real3 getVz() const; + QUATTR real4 to_real4() const; + + QUATTR Quat getConjugate() const; }; - + QUATTR Quat::Quat(){ n = real(1.0); v = real3(); } - - QUATTR Quat::Quat(real4 q){ + + QUATTR Quat::Quat(const real4& q){ n = q.x; v.x = q.y; v.y = q.z; v.z = q.w; } - - QUATTR Quat::Quat(real n, real3 v){ + + QUATTR Quat::Quat(real n, const real3& v){ this->n = n; this->v = v; } - - QUATTR Quat::Quat(real n, real vx, real vy, real vz){ + + QUATTR Quat::Quat(real n, + real vx,real vy,real vz){ this -> n = n; v.x = vx; v.y = vy; - v.z = vz; + v.z = vz; } - - QUATTR Quat Quat::operator+(Quat q){ - return Quat(n+q.n,v+q.v); + + QUATTR Quat Quat::operator+(const Quat& q) const { + return Quat(n+q.n,v+q.v); } - QUATTR void Quat::operator+=(Quat q){ + QUATTR void Quat::operator+=(const Quat& q){ n+=q.n; v+=q.v; } - - QUATTR Quat Quat::operator-(Quat q){ - return Quat(n-q.n,v-q.v); - } - - QUATTR void Quat::operator-=(Quat q){ + + QUATTR Quat Quat::operator-(const Quat& q) const { + return Quat(n-q.n,v-q.v); + } + + QUATTR void Quat::operator-=(const Quat& q){ n-=q.n; v-=q.v; } - - QUATTR Quat Quat::operator*(Quat q){ + + QUATTR Quat Quat::operator*(const Quat& q) const { /* - Product of two quaternions: - q3 = q1*q2 = (n1*n2 - v1*v2, n1*v2 + n2*v1 + v1 x v2) - */ - return Quat(n*q.n-dot(v,q.v),n*q.v+v*q.n+cross(v,q.v)); + Product of two quaternions: + q3 = q1*q2 = (n1*n2 - v1*v2, n1*v2 + n2*v1 + v1 x v2) + */ + return Quat(n*q.n-dot(v,q.v),n*q.v+v*q.n+cross(v,q.v)); } - - QUATTR void Quat::operator*=(Quat q){ - n=n*q.n-dot(v,q.v); - v=n*q.v+v*q.n+cross(q.v,v); + + QUATTR void Quat::operator*=(const Quat& q){ + real n_new = n*q.n-dot(v,q.v); + real3 v_new = n*q.v+v*q.n+cross(v,q.v); + + n = n_new; + v = v_new; } - - QUATTR Quat Quat::operator*(real scalar){ - return Quat(scalar*n,scalar*v); + + QUATTR Quat Quat::operator*(real scalar) const { + return Quat(scalar*n,scalar*v); } - + QUATTR void Quat::operator*=(real scalar){ - n*=scalar; - v*=scalar; + n*=scalar; + v*=scalar; } - - QUATTR Quat operator*(real scalar, Quat q){ - return Quat(scalar*q.n,scalar*q.v); + + QUATTR Quat operator*(real scalar, const Quat& q){ + return Quat(scalar*q.n,scalar*q.v); } - - QUATTR Quat Quat::operator/(real scalar){ - return Quat(n/scalar,v/scalar); + + QUATTR Quat Quat::operator/(real scalar) const { + return Quat(n/scalar,v/scalar); } - + QUATTR void Quat::operator/=(real scalar){ n/=scalar; v/=scalar; } - - QUATTR void Quat::operator=(real4 q){ + + QUATTR void Quat::operator=(const real4& q){ n = q.x; v.x = q.y; v.y = q.z; v.z = q.w; } - - QUATTR real3 Quat::getVx(){ + + QUATTR real3 Quat::getVx() const { real a = n; real b = v.x; real c = v.y; real d = v.z; - real3 vx = make_real3(a*a+b*b-c*c-d*d,2*(b*c+a*d),2*(b*d-a*c)); + real3 vx = make_real3(a*a+b*b-c*c-d*d,real(2.0)*(b*c+a*d),real(2.0)*(b*d-a*c)); return vx; } - - QUATTR real3 Quat::getVy(){ + + QUATTR real3 Quat::getVy() const { real a = n; real b = v.x; real c = v.y; real d = v.z; - return make_real3(2*(b*c-a*d),a*a-b*b+c*c-d*d,2*(c*d+a*b)); + return make_real3(real(2.0)*(b*c-a*d),a*a-b*b+c*c-d*d,real(2.0)*(c*d+a*b)); } - - QUATTR real3 Quat::getVz(){ + + QUATTR real3 Quat::getVz() const { real a = n; real b = v.x; real c = v.y; real d = v.z; - return make_real3(2*(b*d+a*c),2*(c*d-a*b),a*a-b*b-c*c+d*d); + return make_real3(real(2.0)*(b*d+a*c),real(2.0)*(c*d-a*b),a*a-b*b-c*c+d*d); } - - QUATTR real4 Quat::to_real4(){ + + QUATTR real4 Quat::to_real4() const { real4 r4; r4.x = n; r4.y = v.x; @@ -160,32 +164,32 @@ namespace uammd{ r4.w = v.z; return r4; } - + QUATTR real4 make_real4(Quat q){ return q.to_real4(); } - QUATTR Quat Quat::getConjugate(){ + QUATTR Quat Quat::getConjugate() const { return Quat(n,-v); } - /* - Returns the quaternion that encondes a rotation of ang radians - around the axis vrot + /* + Returns the quaternion that encondes a rotation of ang radians + around the axis vrot q = (cos(phi/2),vrot·sin(phi/2)) - */ + */ QUATTR Quat rotVec2Quaternion(real3 vrot, real phi){ real norm = (dot(vrot,vrot)); - if (norm==0) return Quat(1.0,0.,0.,0.); + if (norm==real(0.0)) return Quat(1.0, 0., 0., 0.); vrot*=rsqrt(norm); // The rotation axis must be a unitary vector real cphi2, sphi2; real* cphi2_ptr = &cphi2; real* sphi2_ptr = &sphi2; - sincos(phi*0.5,sphi2_ptr,cphi2_ptr); + sincos(phi*real(0.5),sphi2_ptr,cphi2_ptr); Quat q = Quat(cphi2,sphi2*vrot); return q; } - + QUATTR Quat rotVec2Quaternion(real3 vrot){ // If no angle is given the rotation angle is the modulus of vrot real phi = sqrt(dot(vrot,vrot)); @@ -197,10 +201,10 @@ namespace uammd{ To speed up the computation, we write the rotation in the next form: v' = v + 2*q_n*(q_v x v) + 2*q_v * (q_v x v). See https://fgiesen.wordpress.com/2019/02/09/rotating-a-single-vector-using-a-quaternion/ - */ - QUATTR real3 rotateVector(Quat q, real3 v){ - real3 aux = 2*cross(q.v,v); - return v+q.n*aux+cross(q.v,aux); + */ + QUATTR real3 rotateVector(const Quat& q, const real3& v){ + real3 aux = real(2.0)*cross(q.v,v); + return v+q.n*aux+cross(q.v,aux); } }