Mathematica 9 is now available
Student Support Forum
-----
Student Support Forum: 'need reference for this computation of determin...' topicStudent Support Forum > General > Archives > "need reference for this computation of determin..."

Next Comment >Help | Reply To Topic
Author Comment/Response
Fred
06/18/09 00:19am

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:
It's in FORTRAN, and the construction
sum(F(:,1)*F(:,2))
means
"find the dot product of row 1 with column 2."
(You can swap "row" with "column" because F is symmetric.)


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


URL: ,

Subject (listing for 'need reference for this computation of determin...')
Author Date Posted
need reference for this computation of determin... Fred 06/18/09 00:19am
Re: need reference for this computation of dete... Nooj 08/20/09 5:55pm
Next Comment >Help | Reply To Topic