Project report on
Application of Gaussain quadrature in FEM
Submitted by
Khaja Mohiuddin (2010A2PS353H)
Under supervision of
Dr. Chandu Parimi
Department of Civil Engineering BITS-Pilani Hyderabad Campus November 2012
ACKNOWLEDGEMENT
I would like to express my sincere thanks and gratitude to Dr. Dr. Chandu Parimi for his continuous and unfailing support, guidance and help, which have been invaluable during the course of this project. pro ject. His knowledge, insight and constant motivation at each step of the project has been instrumental in its completion.
S.M.Khaja Mohiuddin (2010A2PS353H)
Table Table of contents 1. 2. 3. 4. 5. 6. 7. 8. 9. 10. 11.
Introducti on Gaussi an quadrature Tri angul ar co-ordinates Shape f unctions Gaussian quadrature f or tri angl es Sampl e code Rectangular elements sample code Tetrahedral elements ele ments Gaussian quadrature for tetrahedral tetrahedral elements sample code
04 05 06 09 11 12 15 17 18 20 21
R ef erences
22
Introduction
The Finite Element Method (FEM) is a numerical technique to find approximate solutions of partial differential equations.It was originated from the need of solving complex elasticity and structural analysis problems in Civil and other engineering branches. FEM allows for detailed visualization and indicates the distribution of stresses and strains inside the body of a structure.This powerful design tool has significantly improved both the standard of engineering designs and the methodology of the design process in many industrial applications. The basic finite element problem can be divided into following steps 1. Discretisation of the continuum: The continuum is divided into a number of elements by imaginary lines or surfaces. 2. Identification of variables: The elements are assumed to be connected at their intersecting points referred to as nodal points. At each node, unknown displacements are to be prescribed. 3. Choice of approximating functions: Displacement function is the starting point of the mathematical mathematical analysis. This represents the variation of the displacement within the element.The displacement function may be approximated in the form a linear function or a higherorder function.A convenient way to express it is by polynomial expressions. The shape or geometry of the element may also be approximated. 4. Formation of the element stiffness matrix: After continuum is d discretised iscretised with desired element shapes, the individual element stiffness matrix is formulated.The geometry of the element is defined in reference to the global frame. Coordinate transformation must be done for elements where it is necessary. nece ssary. 5. Formation of overall stiffness matrix: After the element stiffness matrices in global coordinates are formed, they are assembled to form the overall stiffness matrix. The assembly is done through the nodes nod es which are common to adjacent elements. The overall stiffness matrix is symmetric and banded. 6. Formation of the element loading matrix: The loading forms an essential parameter in any structural engineering problem. prob lem. The loading inside an element is transferred at the nodal points and consistent element matrix is formed. 7. Formation of the overall loading matrix: Like the overall stiffness matrix, the element loading matrices are assembled to form the overall loading matrix. This matrix has one column per loading case and it is either a column vector or a rectangular matrix depending on the number of loading cases. 8. Incorporation of boundary conditions: The boundary restraint conditions are to be imposed in the stiffness matrix. There are various techniques available to satisfy the boundary conditions. 9. Solution Solution of simul simultan taneous eous e uation uations: s: The unknown nodal displacements
are calculated by the multiplication of force vector with the inverse of stiffness matrix. 10. Calculation of stresses or stress-resultants: Nodal displacements are utilized for the calculation of stresses or stress-resultants. The integrations, we generally encounter in finite element methods, are quite complicated and it is not possible to find find a solution to those problems. Exact evaluation of the integral associated to the element matrices and the loading vector is not always possible because of the algebraic complexity of the coefficient of the different different equation (i.e., the stiffness influence coefficients, elasticity matrix, loading functions etc.) etc.) so we use numerical integratio. In the finite element analysis, we face the problem of evaluating the following types of integrations
[k] =
( B T) ( D) ( B) dx dy dz
here B is strain displacement matrix and D is constituen t relationship matrix
to find elemental stiffness matrix [k].Approximat [k ].Approximatee solutions to such problems are possible using certain numerical numerical techniques. Several numerical techniques are available, in mathematics for solving definite integration problems,including, mid-point rule, trapezoidal-rule, Simpson’s 1/3rd rule, Simpson’s 3/8th rule and Gauss Quadrature formula. Among these, Gauss Quadrature technique is most useful one for solving problems in finite element method.
Gauss Quadrature Let us start the discussion with one dimensional gaussain quadrature.The n point gaussain quadrature can be represented as follows for a definite integral from -1 to 1. 1
f ( x) dx = − 1
n
i =
( wi f ( xi) )
1
where xi are quadrature points and w i are weights of those points. The location and weight coefficients of Gauss points are calculated by Legendre polynomials. Hence this method is also sometimes referred as Gauss-Legendre Quadrature method.The summation of these values at n sampling points gives the exact solution of a polynomial integrand of an order up to 2n-1. For example, considering sampling at two Gauss points we can get exact solution for a polynomial of an order (2×2-1) (2×2-1) or 3. The use of more number of Gauss points has no effect on accuracy of results but takes more computation time. The following table gives the values of w i and xi for n = 1 to 6 and the integral should be evaluated from -1 to 1. Numbe umberr of Gauss Gauss Point Point Locat Location ion,, x1 points, n
Weight, w 1
1
0.0
2.0
2
±0.5773502692
1.0
3
0.0
0.8888888889 (=8/9)
±0.7745966692
0.555555 0.5555555556 5556 (=5/9) (=5/9)
±0.3399810436
0.6521451549
±0.861363116
0.3478548451
0.0
0.5688888889
±0.5384693101
0.4786286705
±0.9061798459
0.2369268851
±0.2386191861
0.4679139346
±0.6612093865
0.3607615730
±0.9324695142
0.1713244924
4 5
6
For two dimensional integration problems the above mentioned method can be extended by first evaluating the inner integral, keeping η constant, and then evaluating the outer integral.
1
1
f ( x , y) dx dy = − 1 − 1
n
n
i =
1
j =
( wi w j f ( xi , y j) )
1
similarly for 3-d integration. As we have seen now the integration for a function can be done easily by gaussain quadrature if the integrals are from -1 to 1, but once the integral limits change applying quadrature rules become difficult. We can change the limits for a one dimensional problem directly as follows b
( b − a) f ( x) dx = 2 a
n
i =
( b − a) z i + ( a + b) 2 wi f
1
It is not that easy for a 2-d integral as it may have variable limits.T limits.To simplify the application of quadrature rules to 2-d and 3-d integrals we use natural co-ordinate system.
Natural Co-ordinate system Natural coordinate system is basically a local coordinate system which allows the specification of a point within the element by a set of o f dimensionless numbers whose magnitude never exceeds unity. This coordinate system is found to be very effective in formulating the element properties in finite element formulation. This system is defined in such that the magnitude at nodal points will have unity or zero or a convenient set of fractions. It also facilitates the integration to calculate element stiffness. Triangular elements: The natural coordinate system for a triangular element is generally called as triangular coordinate system. The coordinate of any point P inside the triangle is x,y in Cartesian coordinate system.Here system.Here the co-ordinates ,L1,L2,L3 ,L1,L2,L3 can be used to define the location of points in terms of natural co-ordinates.The point p can be defined by the following set of natural co-ordinates L1 = A1/A ; L2 = A2/A ; L3 = A3/A ; where A1,A2,A3,A are the area of the triangles P23,P13,P12,123 It can be seen that A = A1+A2+A3; so L1+L2+L3 = 1; Therefore, the natural coordinate of three nodes will be: node 1 (1,0,0); node 2 (0,1,0); and node 3 (0,0,1).
Shape functions Shape Function using Area Coordinates Now as we have a understanding about triangular co-ordinate system we can move to shape functions.Shape functions are functions which interpolate the points in the geometry in terms of their nodes.So if N1,N2,N3 N1,N2,N3 are shape functions the we can write x = N1 x1 +N2 x2 +N3 x3 y
= N1 y1 +N2 y2 +N3 y3
Its is important to know that area co-ordinates are dependent and L1+L2+L3 = 1 and also that L1 = 1 while L2,L3=0 at node 1 and L2 = 1 while L1,L3=0 at node 2 and L3 = 1 while L2,L1 =0.So =0.So we can write any point p in vector notations as a linear combination of L1,L2,L3 L1,L2,L3
P = a1L1 + a2L2 + a3L3; and considering nodes x1 x2
x3 thus
=
a1 a2
a3
=
1
0
0
1
0
0
1
0
0
1
0
0
a1 0 a2 1 a3 0
x1 0 x2 1 x3 0
so P = ( L1 L2 L3 )
1
0
0
1
0
0
x1 0 x2 1 x3 0
so we can say see that x = L1 x1 +L2 x2 +L3 x3 y = L1 y1 +L2 y2 +L3 y3 thus the shape functions for a 3 noded triangle are L1,L2,L3 itself . Six noded triangle The below figure shows a triangular element with six nodes. The additional three nodes (4, 5, and 6) are situated at the midpoints of the sides of the element.
Using the above field variable function, one can reach the following expression using interpolation function and the nodal values. x = N1 x1 +N2 x2 +N3 x3 +N4 x4 +N5 x5 +N6 x6 y = N1 y1 +N2 y2 +N3 y3+ N5 y4 +N5 y5 +N6 y6 Here, the every shape function must be such that its value will be unity if evaluated at its related node and zero if evaluated at any of the other five nodes. Now the interpolation functions can be constructed with the help of area coordinates from the above diagram. For example, the interpolation function N should be unity at node 1 and zero at all other five nodes. According to the above diagram, the value of L L1 is 1 at node 1 and ½ at node 4 and 6. Again, L will be 0 at nodes 2, 3 and 5. To satisfy all these conditions, one can propose following expression N1(x,y) = N1(L1,L2,L3) = L1(L1-1/2) Evaluating the above expression, the value of N is becoming ½ at node 1 though it must become unity. unity. Therefore, the above expression is slightly modified satisfying all the conditions and will be as follows: N1(L1,L2,L3) N1(L1,L2,L3) = 2*L1(L1-1/2)=L1(2 L1 - 1) the above equation assures the required conditions at all the six nodes and is a quadratic function, as L is a linear function of x x and y. The remaining five interpolation functions can also be obtained in similar fashion applying the required nodal conditions. Thus, the shape function for the six node triangle element can be written as given below. N1 = L1(2 L1(2 L1 - 1) N2 = L2(2 L2(2 L2 - 1) N3 = L2(2 L2(2 L3 - 1) N4 = 4L1L2 N5 = 4L2L3 N6 = 4L3L1
Gauss Quadrature for Triangular Elements As we have changed the co-ordinate system in triangular elements to triangular co-ordinates the quadrature also changes as follows
= I F ( L1 , L2 , L3) dx dy =
n
i ( wi f ( L1 i =
i
, L2 ,
L3
)) J
i
1
Where,J is jacobian to transform dx,dy and is equal to area for triangular system L terms are the triangular area coordinates and the w terms are the weights associated withthose coordinates. The locations of integration points are shown in figure
For a linear triangle(ie. one sampling point) w1 = 1 L1 = 1/3 L2 = 1/3 L3 = 1/3 For a quadratic qu adratic triangle (ie. n=3) w1 = 1/3 L1 = 1/2 L2 = 1/2 L3 = 0 w1 = 1/3
L1 = 0 L2 = 1/2 L3 = 1/2
w1 = 1/3
L1 = 1/2 L2 = 0 L3 = 1/2
For a cubic triangle (ie.n=7) w1 = 27/60 L1 = 1/3 L2 = 1/3 L3 = 1/3 w1 = 8/60
L1 = 1/2 L2 = 1/2 L3 = 0
w1 = 8/60
L1 = 0 L2 = 1/2 L3 = 1/2
w1 = 8/60
L1 = 1/2 L2 = 0 L3 = 1/2
w1 = 3/60
L1 = 1 L2 = 1 L3 = 0
w1 = 3/60
L1 = 0 L2 = 1 L3 = 1
w1 = 3/60
L1 = 1 L2 = 0 L3 = 1
Code for a 3 node triangle
In the code below x1,y1,x2,y2,x3, x1,y1,x2,y2,x3,y3 y3 n f
are x,y co -ordinates for a 3 noded triangle Number of sampling points Function to be integrated
function out = shape(x1,y1,x2,y2,x3,y3,n,f); if(n==1) if (n==1) mata = [1 1 1;x1 x2 x3;y1 y2 y3]; q = mata*[1/3;1/3;1/3]; mata*[1/3;1/3;1/3]; area = det(mata)/2; det(mata)/2; if(area<0) if (area<0) area = area*(-1); end notarea = f(q(2,1),q(3,1)); out = area*notarea; end if(n==3) if (n==3) mata = [1 1 1;x1 x2 x3;y1 y2 y3]; area = det(mata)/2; det(mata)/2; if(area<0) if (area<0) area = area*(-1); end q1 = mata*[0.5;0.5;0]; q2 = mata*[0;0.5;0.5]; q3 = mata*[0.5;0;0.5]; notarea = (1/3*f(q1(2,1),q1(3,1)))+(1/3*f(q2(2,1),q2(3,1)))+(1/3*f(q3(2,1),q3(3,1) )); out = area*notarea; end if(n==7) if (n==7) mata = [1 1 1;x1 x2 x3;y1 y2 y3]; area = det(mata)/2; det(mata)/2; if(area<0) if (area<0) area = area*(-1); end q1 = mata*[1/3;1/3;1/3]; q2 = mata*[0.5;0.5;0]; q3 = mata*[0;0.5;0.5]; q4 = mata*[0.5;0;0.5]; q5 = mata*[1;0;0]; q6 = mata*[0;1;0]; q7 = mata*[0;0;1]; notarea = (27/60*f(q1(2,1),q1(3,1)))+(8/60*f(q2(2,1),q2(3,1)))+(8/60*f(q3(2,1),q3( 3,1)))+(8/60*f(q4(2,1),q4(3,1)))+(3/60*f(q5(2,1),q5(3,1)))+(3/60*f(q6(2, 1),q6(3,1)))+(3/60*f(q7(2,1),q7(3,1)));
out = area*notarea; end
Code for a 6-noded triangle In the code below x1,y1,x2,y2,x3,y3,x4,y4,x5,y5,x6,y6 are x,y co -ordinates -ordin ates for a 6 noded triangle the order of co-ordinate input is important first 3 points are corner points given in counter clockwise direction and next 3 points are mid points given in counter clockwise direction. n Number of sampling points f Function to be integrated function out = shape6(x1,y1,x2,y2,x3,y3,x4,y4,x5,y5,x6,y6,n,f); if(n==1) if (n==1) mata = [1 1 1 1 1 1;x1 x2 x3 x4 x5 x6;y1 y2 y3 y4 y5 y6]; [n1,n2,n3,n4,n5,n6] = ns(1/3,1/3,1/3); q = mata*[n1;n2;n3;n4;n5;n6] matb = [1 1 1;x1 x2 x3;y1 y2 y3]; area = det(matb)/2; l1=1/3;l2=1/3;l3=1/3; if(area<0) if (area<0) area = area*(-1); end notarea = f(q(2,1),q(3,1)); out = notarea*area; end if(n==3) if (n==3) mata = [1 1 1 1 1 1;x1 x2 x3 x4 x5 x6;y1 y2 y3 y4 y5 y6]; matb = [1 1 1;x1 x2 x3;y1 y2 y3]; area = det(matb)/2; if(area<0) if (area<0) area = area*(-1); end [n1,n2,n3,n4,n5,n6] = ns(0.5,0.5,0); q1 = mata*[n1;n2;n3;n4;n5;n6]; [n1,n2,n3,n4,n5,n6] = ns(0,0.5,0.5); q2 = mata*[n1;n2;n3;n4;n5;n6]; [n1,n2,n3,n4,n5,n6] = ns(0.5,0,0.5); q3 = mata*[n1;n2;n3;n4;n5;n6]; notarea = (1/3*f(q1(2,1),q1(3,1)))+(1/3*f(q2(2,1),q2(3,1)))+(1/3*f(q3(2,1),q3(3,1)) ); out = area*notarea; end if(n==7) if (n==7) mata = [1 1 1 1 1 1;x1 x2 x3 x4 x5 x6;y1 y2 y3 y4 y5 y6]; matb = [1 1 1;x1 x2 x3;y1 y2 y3]; area = det(matb)/2; if(area<0) if (area<0)
area = area*(-1); end [n1,n2,n3,n4,n5,n6] = ns(1/3,1/3,1/3); q1 = mata*[n1;n2;n3;n4;n5;n6]; [n1,n2,n3,n4,n5,n6] = ns(0.5,0.5,0); q2 = mata*[n1;n2;n3;n4;n5;n6]; [n1,n2,n3,n4,n5,n6] = ns(0,0.5,0.5); q3 = mata*[n1;n2;n3;n4;n5;n6]; [n1,n2,n3,n4,n5,n6] = ns(0.5,0,0.5); q4 = mata*[n1;n2;n3;n4;n5;n6]; [n1,n2,n3,n4,n5,n6] = ns(1,0,0); q5 = mata*[n1;n2;n3;n4;n5;n6]; [n1,n2,n3,n4,n5,n6] = ns(0,1,0); q6 = mata*[n1;n2;n3;n4;n5;n6]; [n1,n2,n3,n4,n5,n6] = ns(0,0,1); q7 = mata*[n1;n2;n3;n4;n5;n6]; notarea = (27/60*f(q1(2,1),q1(3,1)))+(8/60*f(q2(2,1),q2(3,1)))+(8/60*f(q3(2,1),q3(3 ,1)))+(8/60*f(q4(2,1),q4(3,1)))+(3/60*f(q5(2,1),q5(3,1)))+(3/60*f(q6(2,1) ,q6(3,1)))+(3/60*f(q7(2,1),q7(3,1))); out = area*notarea; end end function [n1,n2,n3,n4,n5,n6]= ns(l1,l2,l3) n1 = l1*(2*l1-1); n2 = l2*(2*l2-1); n3 = l3*(2*l3-1); n4 = 4*l1*l2; n5 = 4*l2*l3; n6 = 4*l3*l1; end
Rectangular Rectangular elements It is easy to define the new n ew co-ordinate system for rectangular elements. Just a change of origin will do the job. j ob. Let the co-ordinates in new system sys tem be η and ξ.
Now we define the interpolation function of the element in natural coordinate system and can be found as:
N1 N2 = Ni := N3 N4
(1 − (1 + (1 +
(1
−
) ) ) )
(1
)
− 4
(1
)
− 4
(1
)
+ 4
(1
)
+ 4
These can found by using simple lagrange interpolation theorem
So now x and y can be written as x = N1 x1 x1 + N2 x2
+
N3 x3
+
N4 x4
y = N1 y1 y1 + N2 y2
+
N3 y3
+
N4 y4
We have got x,y x,y in terms of nodes now we should convert dxdy to dηdξ this can be done by the introduction of jacobian term defined as follows
So I can be evaluated using gaussain quadrature as follows
1
n
1
I= − 1 − 1
(
,
) J d d
=
n
i =
1
j =
1
( wi w j)
(
,
) J
Sample code for a four fo ur noded rectangular rectangu lar element element Here x1,y1,x2,y2,x3,y3,x4,y4 are co ordinates ordinate s of the given rectangular element in anti-clockwise direction n Number of points f Function to be integrated function out = shapequad(x1,y1,x2,y2,x3,y3,x4,y4,n,f); [r,cnst] = lgwt(n,-1,1); integretion = 0; for i=1:n jx=0; l1= r(i,1); for j=1:n l2 = r(j,1); [n1,n2,n3,n4,jaco]=ns(l1,l2,x1,y1,x2,y2,x3,y3,x4,y4); x=n1*x1+n2*x2+n3*x3+n4*x4; y=n1*y1+n2*y2+n3*y3+n4*y4; q = jaco*f(x,y); jx = jx + (cnst(j,1)*q); end integretion = integretion + (cnst(i,1)*jx); end out = integretion; integretion; end function [n1,n2,n3,n4,j] = ns(l1,l2,x1,y1,x2,y2,x3,y3,x4,y4) n1 = 0.25*(1-l1)*(1-l2); n2 = 0.25*(1+l1)*(1-l2); n3 = 0.25*(1+l1)*(1+l2); n4 = 0.25*(1-l1)*(1+l2); op=[-0.25*(1-l1) 0.25*(1-l1) 0.25*(1+l1) -0.25*(1+l1);-0.25*(1-l2) -0.25*(1+l2) 0.25*(1+l2) 0.25*(1-l2)]; oq = [x1 y1;x2 y2;x3 y3;x4 y4]; j = op*oq; j = det(j); end
Tetrahedral etrahed ral elements elemen ts There are two basic families of three-dimensional elements similar to two dimensional case.Extension case.Extension of triangular elements will produce tetrahedrons in three dimensions. Similarly, Similarly, rectangular parallelepipeds are generated on the extension of rectangular elements. elements. Figure shows few commonly used solid elements for finite element analysis.
Tetrahedral elements The simplest element of the tetrahedral family is a four node tetrahedron as shown in figure. The node numbering has been followed in sequential manner, manner, i.e, in this case anti-clockwise direction. Similar to the area coordinates, the concept of volume coordinates has been introduced here. The coordinates of the x, y, and z) nodes are defined both in Cartesian and volume coordinates. Point P( x, as shown in below figure is an arbitrary point in the tetrahedron.
The linear shape function for this element can be expressed as, NT = [ L1 L2 L3 L4] Here L1,L2,L3,L4 L1,L2,L3,L4 are natural co ordinates defined as follows
Li
=
Vi V
Where Vi is the volume of the sub element which is bound by point P and face i and V is the total volume of o f the element. For example L i may be interpreted as the ratio of the volume of the sub element P234 to the total volume of the element 1234. The volume of the element V is given by the determinant of the nodal coordinates as follows:
V
:=
1 1 x1 6 y1 z1
1
1
x2 x3 y2 y3 z2 z3
x4 y4 z4 1
The relationship between the Cartesian and natural coordinates of point P may be expressed as
1 L1 1 x := L2 x1 y L3 y1 z L4 z1
1
1
x2 x3 y2 y3 z2 z3
x4 y4 z4 1
so L1 L2 L3 L4 can found by inversing the matrix
Gaussain quadrature for tetrahedral elements The Gauss Quadrature for triangles can be effectively extended to include
tetrahedron elements in terms of tetrahedron volume coordinates.
I = F ( L1 , L2 , L3 , L4) dx dy dz =
n
i ( wi f ( L1 i =
,
i
i
L2 , L3 , L4
1
J
Where, L terms are the volume coordinates and the wi terms are the weights associated with those coordinates. The locations locations of Gauss points are shown in figure below
The sampling points and their associated weights are described below: For sampling point = 1 (Linear tetrahedron) w1 = 1
L11=L21=L31=L41=1/4
For sampling point = 4 (Quadratic tetrahedron) w1 = 1/4
L11=0.5854102 L21=L31=L41=0.1381966
w2 = 1/4
L22=0.5854102 L12=L32=L42=0.1381966
w3 = 1/4
L33=0.5854102 L13=L23=L43=0.1381966
w4 = 1/4
L44=0.5854102 L24=L34=L14=0.1381966
For sampling points = 5 (Cubic tetrahedron) w1 = -4/5
L11=L21=L31=L41=1/4
w2 = 9/20
L12=1/3 L22=L32=L42=1/6
w3 = 9/20
L24=1/3 L14=L34=L44=1/6
w4 = 9/20
L33=1/3 L13=L23=L43=1/6
w5 = 9/20
L44=1/3 L24=L34=L14=1/6
Sample code for a tetrahedron
))
i
function out = shapetetra(x1,y1,z1,x2,y2,z2,x3,y3,z3,x4,y4,z4,n,f); if(f==1) if (f==1) mata = [1 1 1 1;x1 x2 x3 x4;y1 y2 y3 y4;z1 z2 z3 z4]; volume = det(mata)/6; det(mata)/6; if(volume<0) if (volume<0) volume = volume*(-1); volume*(-1); end out = volume; else if(n==1) if (n==1) mata = [1 1 1 1;x1 x2 x3 x4;y1 y2 y3 y4;z1 z2 z3 z4]; q = mata*[1/4;1/4;1/4;1/4]; volume = det(mata)/6; if(volume<0) if (volume<0) volume = volume*(-1); volume*(-1); end notvolume = f(q(2,1),q(3,1),q(4,1)); out = volume*notvolume; end if(n==4) if (n==4) mata = [1 1 1 1;x1 x2 x3 x4;y1 y2 y3 y4;z1 z2 z3 z4]; volume = det(mata)/6; if(volume<0) if (volume<0) volume = volume*(-1); volume*(-1); end q1 = mata*[0.5854102;0.1381966;0.1381966;0.1381966]; q2 = mata*[0.1381966;0.5854102;0.1381966;0.1381966]; q3 = mata*[0.1381966;0.1381966;0.5854102;0.1381966]; q4 = mata*[0.1381966;0.1381966;0.1381966;0.5854102]; notvolume = (1/4*f(q1(2,1),q1(3,1),q1(4,1)))+(1/4*f(q2(2,1),q2(3,1),q2(4,1)))+(1/4*f (q3(2,1),q3(3,1),q3(4,1)))+(1/4*f(q4(2,1),q4(3,1),q4(4,1))); out = volume*notvolume; end if(n==5) if (n==5) mata = [1 1 1 1;x1 x2 x3 x4;y1 y2 y3 y4;z1 z2 z3 z4]; volume = det(mata)/6; if(volume<0) if (volume<0) volume = volume*(-1); volume*(-1); end q1 = mata*[1/4;1/4;1/4;1/4]; q2 = mata*[1/3;1/6;1/6;1/6]; q3 = mata*[1/6;1/3;1/6;1/6]; q4 = mata*[1/6;1/6;1/3;1/6]; q5 = mata*[1/6;1/6;1/6;1/3]; notvolume = ((-4/5)*f(q1(2,1),q1(3,1),q1(4,1)))+((9/20)*f(q2(2,1),q2(3,1),q2(4,1)))+ ((9/20)*f(q3(2,1),q3(3,1),q3(4,1)))+(9/20*f(q4(2,1),q4(3,1),q4(4,1)))+(9 /20*f(q5(2,1),q5(3,1),q5(4,1))); out = volume*notvolume; end end
References: Wikipedia
Google Nptel video lectures Introduction to Finite Element Methods course at University of colorado at Boulder
