Skip to content
GitLab
Menu
Projects
Groups
Snippets
Help
Help
Support
Community forum
Keyboard shortcuts
?
Submit feedback
Contribute to GitLab
Sign in / Register
Toggle navigation
Menu
Open sidebar
Development
openfoam
Commits
661d1f53
Commit
661d1f53
authored
Apr 25, 2019
by
Mark Olesen
Committed by
Andrew Heather
Apr 25, 2019
Browse files
ENH: additional rotation tests (
#1292
)
parent
bc53e50c
Changes
2
Hide whitespace changes
Inline
Side-by-side
applications/test/quaternion/Make/options
View file @
661d1f53
/* EXE_INC = */
/* EXE_LIBS = */
EXE_INC = \
-I$(LIB_SRC)/meshTools/lnInclude
EXE_LIBS = \
-lmeshTools
applications/test/quaternion/Test-quaternion.C
View file @
661d1f53
...
...
@@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd |
\\ / A nd |
Copyright (C) 2019 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
| Copyright (C) 2011-2016 OpenFOAM Foundation
...
...
@@ -34,12 +34,29 @@ Description
#include
"argList.H"
#include
"quaternion.H"
#include
"septernion.H"
#include
"
mathematicalConstants
.H"
#include
"
unitConversion
.H"
#include
"Tuple2.H"
#include
"IOstreams.H"
#include
"transform.H"
#include
"axisAngleRotation.H"
#include
"EulerCoordinateRotation.H"
using
namespace
Foam
;
using
namespace
Foam
::
coordinateRotations
;
void
printRotation
(
const
tensor
&
rot
)
{
Info
<<
"rotation = "
<<
rot
<<
nl
;
}
void
printRotation
(
const
quaternion
&
quat
)
{
printRotation
(
quat
.
R
());
Info
<<
"quaternion "
<<
quat
<<
nl
;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
...
...
@@ -57,7 +74,24 @@ int main(int argc, char *argv[])
"vector"
,
"Rotate by '(yaw pitch roll)' in degrees"
);
argList
::
addOption
(
"euler"
,
"vector"
,
"Rotate by '(phi theta psi)' in degrees"
);
argList
::
addOption
(
"xyz"
,
"vector"
,
"Rotate about x-y-z axes in degrees"
);
argList
::
addOption
(
"zyx"
,
"vector"
,
"Rotate about z-y-x axes in degrees"
);
argList
::
addOption
(
"rotate-angle"
,
...
...
@@ -70,34 +104,77 @@ int main(int argc, char *argv[])
if
(
args
.
readIfPresent
(
"rollPitchYaw"
,
rotVector
))
{
Info
<<
"Rotate by"
<<
nl
Info
<<
nl
<<
"Rotate by"
<<
nl
<<
" roll "
<<
rotVector
.
x
()
<<
nl
<<
" pitch "
<<
rotVector
.
y
()
<<
nl
<<
" yaw "
<<
rotVector
.
z
()
<<
nl
;
// degToRad
rotVector
*=
constant
::
mathematical
::
pi
/
180
.
0
;
rotVector
*=
degToRad
();
const
quaternion
quat
(
quaternion
::
rotationSequence
::
XYZ
,
rotVector
);
Info
<<
"quaternion "
<<
quat
<<
endl
;
Info
<<
"rotation = "
<<
quat
.
R
()
<<
endl
;
printRotation
(
quat
);
}
if
(
args
.
readIfPresent
(
"yawPitchRoll"
,
rotVector
))
{
Info
<<
"Rotate by"
<<
nl
Info
<<
nl
<<
"Rotate by"
<<
nl
<<
" yaw "
<<
rotVector
.
x
()
<<
nl
<<
" pitch "
<<
rotVector
.
y
()
<<
nl
<<
" roll "
<<
rotVector
.
z
()
<<
nl
;
// degToRad
rotVector
*=
constant
::
mathematical
::
pi
/
180
.
0
;
rotVector
*=
degToRad
();
const
quaternion
quat
(
quaternion
::
rotationSequence
::
ZYX
,
rotVector
);
Info
<<
"quaternion "
<<
quat
<<
endl
;
Info
<<
"rotation = "
<<
quat
.
R
()
<<
endl
;
printRotation
(
quat
);
}
if
(
args
.
readIfPresent
(
"euler"
,
rotVector
))
{
Info
<<
nl
<<
"Rotate by"
<<
nl
<<
" phi "
<<
rotVector
.
x
()
<<
nl
<<
" theta "
<<
rotVector
.
y
()
<<
nl
<<
" psi "
<<
rotVector
.
z
()
<<
nl
;
printRotation
(
euler
(
rotVector
,
true
).
R
());
}
if
(
args
.
readIfPresent
(
"xyz"
,
rotVector
))
{
Info
<<
nl
<<
"Rotate about"
<<
nl
<<
" x "
<<
rotVector
.
x
()
<<
nl
<<
" y "
<<
rotVector
.
y
()
<<
nl
<<
" z "
<<
rotVector
.
z
()
<<
nl
;
Vector
<
tensor
>
Rs
(
axisAngle
(
vector
(
1
,
0
,
0
),
rotVector
.
x
(),
true
).
R
(),
axisAngle
(
vector
(
0
,
1
,
0
),
rotVector
.
y
(),
true
).
R
(),
axisAngle
(
vector
(
0
,
0
,
1
),
rotVector
.
z
(),
true
).
R
()
);
printRotation
(
Rs
.
x
()
&
Rs
.
y
()
&
Rs
.
z
());
}
if
(
args
.
readIfPresent
(
"zyx"
,
rotVector
))
{
Info
<<
nl
<<
"Rotate about"
<<
nl
<<
" z "
<<
rotVector
.
x
()
<<
nl
<<
" y "
<<
rotVector
.
y
()
<<
nl
<<
" x "
<<
rotVector
.
z
()
<<
nl
;
Vector
<
tensor
>
Rs
(
axisAngle
(
vector
(
1
,
0
,
0
),
rotVector
.
z
(),
true
).
R
(),
axisAngle
(
vector
(
0
,
1
,
0
),
rotVector
.
y
(),
true
).
R
(),
axisAngle
(
vector
(
0
,
0
,
1
),
rotVector
.
x
(),
true
).
R
()
);
printRotation
(
Rs
.
z
()
&
Rs
.
y
()
&
Rs
.
x
());
}
if
(
args
.
found
(
"rotate-angle"
))
{
const
Tuple2
<
vector
,
scalar
>
axisAngle
...
...
@@ -105,60 +182,61 @@ int main(int argc, char *argv[])
args
.
lookup
(
"rotate-angle"
)()
);
Info
<<
"Rotate"
<<
nl
Info
<<
nl
<<
"Rotate"
<<
nl
<<
" about "
<<
axisAngle
.
first
()
<<
nl
<<
" angle "
<<
axisAngle
.
second
()
<<
nl
;
const
quaternion
quat
(
axisAngle
.
first
(),
axisAngle
.
second
()
*
constant
::
mathematical
::
pi
/
180
.
0
// degToRad
degToRad
(
axisAngle
.
second
()
)
);
Info
<<
"quaternion "
<<
quat
<<
endl
;
Info
<<
"rotation = "
<<
quat
.
R
()
<<
endl
;
printRotation
(
quat
);
Info
<<
"transform Ra = "
<<
Ra
(
axisAngle
.
first
()
/
mag
(
axisAngle
.
first
()),
axisAngle
.
second
()
*
constant
::
mathematical
::
pi
/
180
.
0
)
<<
end
l
;
)
<<
n
l
;
Info
<<
"-ve Ra = "
<<
Ra
(
axisAngle
.
first
()
/
mag
(
axisAngle
.
first
()),
-
axisAngle
.
second
()
*
constant
::
mathematical
::
pi
/
180
.
0
)
<<
end
l
;
)
<<
n
l
;
}
Info
<<
nl
<<
nl
;
quaternion
q
(
vector
(
1
,
2
,
3
),
0
.
7853981
);
Info
<<
"q "
<<
q
<<
end
l
;
Info
<<
"q "
<<
q
<<
n
l
;
vector
v
(
0
.
1
,
0
.
4
,
2
.
1
);
Info
<<
"v "
<<
v
<<
end
l
;
Info
<<
"v "
<<
v
<<
n
l
;
Info
<<
"inv(q)*q "
<<
inv
(
q
)
*
q
<<
end
l
;
Info
<<
"inv(q)*q "
<<
inv
(
q
)
*
q
<<
n
l
;
Info
<<
"q*quaternion(0, v)*conjugate(q) "
<<
q
*
quaternion
(
0
,
v
)
*
conjugate
(
q
)
<<
end
l
;
<<
q
*
quaternion
(
0
,
v
)
*
conjugate
(
q
)
<<
n
l
;
Info
<<
"q.R() "
<<
q
.
R
()
<<
end
l
;
Info
<<
"q.transform(v) "
<<
q
.
transform
(
v
)
<<
end
l
;
Info
<<
"q.R() & v "
<<
(
q
.
R
()
&
v
)
<<
end
l
;
Info
<<
"q.R() "
<<
q
.
R
()
<<
n
l
;
Info
<<
"q.transform(v) "
<<
q
.
transform
(
v
)
<<
n
l
;
Info
<<
"q.R() & v "
<<
(
q
.
R
()
&
v
)
<<
n
l
;
Info
<<
"quaternion(q.R()).transform(v) "
<<
(
quaternion
(
q
.
R
()).
transform
(
v
))
<<
end
l
;
<<
(
quaternion
(
q
.
R
()).
transform
(
v
))
<<
n
l
;
Info
<<
"q.invTransform(v) "
<<
q
.
invTransform
(
v
)
<<
end
l
;
Info
<<
"q.invTransform(v) "
<<
q
.
invTransform
(
v
)
<<
n
l
;
septernion
tr
(
vector
(
0
,
0
.
1
,
0
),
q
);
Info
<<
"tr "
<<
tr
<<
end
l
;
Info
<<
"tr "
<<
tr
<<
n
l
;
Info
<<
"inv(tr)*tr "
<<
inv
(
tr
)
*
tr
<<
end
l
;
Info
<<
"inv(tr)*tr "
<<
inv
(
tr
)
*
tr
<<
n
l
;
Info
<<
"tr.transform(v) "
<<
tr
.
transformPoint
(
v
)
<<
end
l
;
Info
<<
"tr.transform(v) "
<<
tr
.
transformPoint
(
v
)
<<
n
l
;
vector
origin
(
1
,
2
,
4
);
...
...
@@ -168,9 +246,9 @@ int main(int argc, char *argv[])
<<
" "
<<
septernion
(
-
origin
)
.
transformPoint
(
q
.
transform
(
septernion
(
origin
).
transformPoint
(
v
)))
<<
end
l
;
<<
n
l
;
Info
<<
"Test conversion from and to Euler-angles"
<<
end
l
;
Info
<<
"Test conversion from and to Euler-angles"
<<
n
l
;
vector
angles
(
0
.
1
,
0
.
2
,
0
.
3
);
for
(
int
rs
=
quaternion
::
ZYX
;
rs
<
quaternion
::
XZX
;
++
rs
)
{
...
...
@@ -188,7 +266,7 @@ int main(int argc, char *argv[])
FatalErrorInFunction
<<
"Inconsistent conversion for rotation sequence "
<<
rs
<<
exit
(
FatalError
)
<<
end
l
;
<<
n
l
;
}
}
...
...
@@ -202,9 +280,9 @@ int main(int argc, char *argv[])
ss
[
2
]
=
septernion
(
vector
(
0
,
0
.
3
,
0
),
quaternion
(
0
.
3
,
vector
(
3
,
2
,
1
)));
w
[
2
]
=
0
.
4
;
Info
<<
"average(ss, w) "
<<
average
(
ss
,
w
)
<<
end
l
;
Info
<<
"average(ss, w) "
<<
average
(
ss
,
w
)
<<
n
l
;
Info
<<
"End
\n
"
<<
endl
;
Info
<<
"
\n
End
\n
"
<<
endl
;
return
0
;
}
...
...
Write
Preview
Supports
Markdown
0%
Try again
or
attach a new file
.
Attach a file
Cancel
You are about to add
0
people
to the discussion. Proceed with caution.
Finish editing this message first!
Cancel
Please
register
or
sign in
to comment