Programming Finite Element Method

Question

I am trying to teach myself how about Finite Element Methods.

All of my code is adapted from the following link pages 16-20 http://homepages.cae.wisc.edu/~suresh/ME964Website/M964Notes/Notes/introfem.pdf

I am programming along in Matlab to perform a finite element analysis on a single 8 node cube element. I have defined the xi,eta,zeta local axes (we can think about this as x, y, z for now), so I get the following shape functions:

%%shape functions

zeta = 0:.01:1;
eta = 0:.01:1;
xi = 0:.01:1;

N1 = 1/8*(1-xi).*(1-eta).*(1-zeta);
N2 = 1/8*(1+xi).*(1-eta).*(1-zeta);
N3 = 1/8*(1+xi).*(1+eta).*(1-zeta);
N4 = 1/8*(1-xi).*(1+eta).*(1-zeta);
N5 = 1/8*(1-xi).*(1-eta).*(1+zeta);
N6 = 1/8*(1+xi).*(1-eta).*(1+zeta);
N7 = 1/8*(1+xi).*(1+eta).*(1+zeta);
N8 = 1/8*(1-xi).*(1+eta).*(1+zeta);

The [N] Matrix is to be arranged like this according to the text I am reading:

%N Matrix
N= [N1 0 0 N2 0 0 N3 0 0 N4 0 0 N5 0 0 N6 0 0 N7 0 0 N8 0 0;
    0 N1 0 0 N2 0 0 N3 0 0 N4 0 0 N5 0 0 N6 0 0 N7 0 0 N8 0;
    0 0 N1 0 0 N2 0 0 N3 0 0 N4 0 0 N5 0 0 N6 0 0 N7 0 0 N8];

To find the [B] matrix i have to use the following [D] matrix:

%%Del Matrix for node i
%[  d/dx   0     0
%    0    d/dy   0
%    0     0    d/dz        . . .
%   d/dy  d/dx   0
%    0    d/dz  d/dy
%   d/dz   0    d/dx   ]

which is an operator to go on [N]. (B=DN)

Later on, as the text shows, I will be making calculations involving integrals of this [B] matrix over the volume of this element.

So, my question is, how can I store these polynomial shape functions in a matrix, operate on them with differentiation, and then integrate them numerically. I can tell with the way I have this set up right now, that it wont work because I have defined the functions as a vector over an interval [0,1] and then storing these vectors in the [N] matrix. Then using diff() function to differentiate appropriately to find the [B] matrix. But since the matrix elements of [B] are now vectors over an interval [0,1] I think that is going to cause problems. How would you guys go about these calculations described in the textbook I posted above?

Answer 1

Solved my problem using anonymous functions and storing the polynomials in a symbolic matrix. example:

syms xi eta zeta
N1= ... %type out in terms of xi eta and zeta
.
.
.
dN1dXi = diff(N1,xi) %symbolic differentiation with respect to xi

can also perform symbolic integration when needed:

intN1 = int(N1,xi,lowerLimit,upperLimit) %symbolic integration with respect to xi

and when ready to substitute in actual values to evaluate the symbolic functions:

subs(N1,{xi,eta,zeta},{value1,value2,value3})

Answer 2

You should check page 24 about how to map from a parametric domain ([0,1]^) to the physical domain.