**A**=<

**pp**> where

**p**is the orientation of a single fiber as shown in the figure below, and the angular brackets denote an average over the distribution function. Although, I explicitly refer to the fiber orientation tensor in this blog, the discussion below is valid for any symmetric tensor, such as a stress tensor, for example.

For non-simple cases, it can be tricky to mentally visualize the tensor. However, if we calculate the eigenvalues and eigenvectors of the matrix

**A**, visualizing the distribution can be done easily by drawing an ellipsoid. The ellipsoid points in a direction determined by the eigenvectors and the shape and size of the ellipsoid is determined by the eigenvalues.

I modified some code (

**VizOrient.m**) that I found online. It works fine with any relatively recent version of Matlab and Octave 3.0 and higher. Below are five examples:

**A=[0.45 0 0; 0 0.45 0; 0 0 0.1];**

VizOrient(A);

VizOrient(A);

**A=[0.90 0 0; 0 0.05 0; 0 0 0.05];**

VizOrient(A);

VizOrient(A);

**A=[0.33 0 0; 0 0.33 0; 0 0 0.33];**

VizOrient(A);

VizOrient(A);

**A=[0.33 0.1 0; 0.1 0.33 0; 0 0 0.33];**

VizOrient(A);

VizOrient(A);

**A=[0.45 0.2 0; 0.2 0.45 0; 0 0 0.1];**

VizOrient(A);

VizOrient(A);

The only problem with Octave is actually a problem with gnuplot. The statement "axis equal" is not strictly implemented, and hence some mental gymnastics could be required. Works great on Matlab, though.

References:

1. Chung and Kwon, Korea-Australia Rheology Journal, 14(4), 2002, 175-188

Code:

**function VizOrient(A)**

%

% take a 3x3 fiber orientation (symmetric) matrix A, and visualize the

% resulting ellipsoid. using a program from Nima Moshtagh, and erasing the 2D

% lines - also correcting the meaning of a, b, and c axes (1/a)

%

%------------------------------------------------------------------------------

% Copyright (c) 2009, Nima Moshtagh

% All rights reserved.

%

% Redistribution and use in source and binary forms, with or without

% modification, are permitted provided that the following conditions are

% met:

%

% * Redistributions of source code must retain the above copyright

% notice, this list of conditions and the following disclaimer.

% * Redistributions in binary form must reproduce the above copyright

% notice, this list of conditions and the following disclaimer in

% the documentation and/or other materials provided with the distribution

%

%THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"

%AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE

%IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE

%ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE

%LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR

%CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF

%SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS

%INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN

%CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)

%ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE

%POSSIBILITY OF SUCH DAMAGE.

%------------------------------------------------------------------------------

N = 20;

[U D V] = svd(A);

%----------------------------------

% generate the ellipsoid at (0,0,0)

%----------------------------------

a = D(1,1);

b = D(2,2);

c = D(3,3);

[X,Y,Z] = ellipsoid(0,0,0,a,b,c,N);

%-------------------------

% rotate and the ellipsoid

%-------------------------

XX = zeros(N+1,N+1);

YY = zeros(N+1,N+1);

ZZ = zeros(N+1,N+1);

for k = 1:length(X),

for j = 1:length(X),

point = [X(k,j) Y(k,j) Z(k,j)]';

P = V * point;

XX(k,j) = P(1);

YY(k,j) = P(2);

ZZ(k,j) = P(3);

end

end

%-----------------

% Plot the ellipse

%-----------------

surf(XX,YY,ZZ);

xlabel('x');

ylabel('y');

zlabel('z');

axis equal

hidden on

end

%

% take a 3x3 fiber orientation (symmetric) matrix A, and visualize the

% resulting ellipsoid. using a program from Nima Moshtagh, and erasing the 2D

% lines - also correcting the meaning of a, b, and c axes (1/a)

%

%------------------------------------------------------------------------------

% Copyright (c) 2009, Nima Moshtagh

% All rights reserved.

%

% Redistribution and use in source and binary forms, with or without

% modification, are permitted provided that the following conditions are

% met:

%

% * Redistributions of source code must retain the above copyright

% notice, this list of conditions and the following disclaimer.

% * Redistributions in binary form must reproduce the above copyright

% notice, this list of conditions and the following disclaimer in

% the documentation and/or other materials provided with the distribution

%

%THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"

%AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE

%IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE

%ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE

%LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR

%CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF

%SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS

%INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN

%CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)

%ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE

%POSSIBILITY OF SUCH DAMAGE.

%------------------------------------------------------------------------------

N = 20;

[U D V] = svd(A);

%----------------------------------

% generate the ellipsoid at (0,0,0)

%----------------------------------

a = D(1,1);

b = D(2,2);

c = D(3,3);

[X,Y,Z] = ellipsoid(0,0,0,a,b,c,N);

%-------------------------

% rotate and the ellipsoid

%-------------------------

XX = zeros(N+1,N+1);

YY = zeros(N+1,N+1);

ZZ = zeros(N+1,N+1);

for k = 1:length(X),

for j = 1:length(X),

point = [X(k,j) Y(k,j) Z(k,j)]';

P = V * point;

XX(k,j) = P(1);

YY(k,j) = P(2);

ZZ(k,j) = P(3);

end

end

%-----------------

% Plot the ellipse

%-----------------

surf(XX,YY,ZZ);

xlabel('x');

ylabel('y');

zlabel('z');

axis equal

hidden on

end

## No comments:

Post a Comment