can anyone recognize and give the name of the algorithm used in this snippet?
i posted this to some math forums, but got no response. perhaps one of you knows the answer. it is not a FORTRAN
question.
SHORT VERSION:
I've got some code that finds a determinant of a matrix F by
taking the
square root of an expression of the form (e*g-f*f), where e,g, and f are
specified elements of
(F^T)*F.
I've never seen this before. Have you? Can you provide me a reference?
LONG
VERSION:
Suppose you have an elastic cube. Suppose you gently deform it by pushing
and pulling on the
corners. Just bend it a little bit; don't jump on it
until it pops or anything like that.
Choose one
of the sides, and a point on that side.
Use "X" to denote the original position of the point (before you
deformed
it), and "x" to denote the new position.
Suppose you determine the deformation gradient F of
the face mapping that
takes X into x.
If that last sentence (or anything else) didn't make sense to you,
ignore
all of the above and suppose you have a symmetric, positive definite matrix
F. If you're still
confused, suppose you have a nice 3x3 matrix F.
Suppose you want the determinant of F. What would you
do? This is what
someone else did:
if ( (face.eq.BOTTOM_FACE)
.or.(face.eq.TOP_FACE))
then
ee = sum(F(:,1)*F(:,1))
ff = sum(F(:,1)*F(:,2))
gg =
sum(F(:,2)*F(:,2))
else if ((face.eq.INTERIOR_FACE)
.or.(face.eq.EXTERIOR_FACE))
then
ee =
sum(F(:,1)*F(:,1))
ff = sum(F(:,1)*F(:,3))
gg = sum(F(:,3)*F(:,3))
else if
((face.eq.RIGHT_FACE)
.or.(face.eq.LEFT_FACE))
then
ee = sum(F(:,2)*F(:,2))
ff =
sum(F(:,2)*F(:,3))
gg = sum(F(:,3)*F(:,3))
else
write(*,*) "ERROR: face not found"
endif
! Jacobian of face mapping
detJb = sqrt(ee*gg-ff*ff)
Anyone know what formula this is?
Can you provide a reference?
He then goes on to use this to find the normal:
if
(face.eq.BOTTOM_FACE) then
nor(1) = F(2,2)*F(3,1) - F(3,2)*F(2,1)
nor(2) = F(3,2)*F(1,1) - F(1,2)*F(3,1)
nor(3) = F(1,2)*F(2,1) - F(2,2)*F(1,1)
else if(face.eq.INTERIOR_FACE) then
nor(1) = F(2,1)*F(3,3) -
F(3,1)*F(2,3)
nor(2) = F(3,1)*F(1,3) - F(1,1)*F(3,3)
nor(3) = F(1,1)*F(2,3) - F(2,1)*F(1,3)
else
if(face.eq.RIGHT_FACE) then
nor(1) = F(2,2)*F(3,3) - F(3,2)*F(2,3)
nor(2) = F(3,2)*F(1,3) - F(1,2)*F(3,3)
nor(3) = F(1,2)*F(2,3) - F(2,2)*F(1,3)
else if(face.eq.EXTERIOR_FACE) then
nor(1) = -F(2,1)*F(3,3) +
F(3,1)*F(2,3)
nor(2) = -F(3,1)*F(1,3) + F(1,1)*F(3,3)
nor(3) = -F(1,1)*F(2,3) + F(2,1)*F(1,3)
else
if(face.eq.LEFT_FACE) then
nor(1) = -F(2,2)*F(3,3) + F(3,2)*F(2,3)
nor(2) = -F(3,2)*F(1,3) +
F(1,2)*F(3,3)
nor(3) = -F(1,2)*F(2,3) + F(2,2)*F(1,3)
else if(face.eq.TOP_FACE) then
nor(1) =
-F(2,2)*F(3,1) + F(3,2)*F(2,1)
nor(2) = -F(3,2)*F(1,1) + F(1,2)*F(3,1)
nor(3) = -F(1,2)*F(2,1) +
F(2,2)*F(1,1)
endif
! make the normal a unit vector
tmp = sqrt( sum(nor(:)*nor(:)) )
nor(:) = nor(:)/tmp
Reference?
Thanks bunches!
Nooj
| |
Re: can anyone recognize and give the name of the algorithm used in this snippet?
I don't think it calculates the determinant; after all, it does not even
take into account all entries of the matrix.
I believe it does this: Let's devide the matrix into the three
vectors a = F(:,1), b=F(:,2) and c=F(:,3). Depending on the variable "face", it calculates the length of the cross
product for a pair of these vectors. For example, if face.eq.BOTTOM_FACE, the norm |a x b| is calculated (where x is the
cross product). Just observe that you can rewrite (|a x b|)^2 as a^2*b^2-(a*b)^2, where a*b is the dot product of a and
b.
It then calculates the vector a x b itself: nor = a x b. This, of course, makes the first part of the code
somewhat redundant as you could as well just calculate the norm of nor.
edit: typo
| |
Re: can anyone recognize and give the name of the algorithm used in this snippet?
tom_p: What you said is correct; this code is known to be accurate and inefficient.
I have found a reference:
It can be shown to be a variant of Nanson's formula ( nda = JF^(-T)NdA ),
and is found in the section "Area of a Curved Surface" on pp 142--144 of
M.D. Greenberg's _Foundations of Applied Mathematics_, Prentice-Hall, 1978.
See in particular formulas (8.45a), (8.45b), and (8.48) on p 143.
Thanks to everyone who pondered this with me!
- Nooj
| | |
