can anyone recognize and give the name of the algorithm used in this snippet?

[nooj]

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

[tom_p]

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

[nooj]

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