I want to evaluate the Surface force integral in an FEM procedure.enter image description here

The basic reference tet is shown in the figure. The faces are numbered corresponding to the node opposite to them. For example the nodes 2-3-4 which forms a face opposite to node 1 is labeled as 1, and so on... I now transform the tet into the volume coordinate denoted by $(\xi_1, \xi_2, \xi_3, \xi_4 )$ which are also the shape functions of the linear element. The surface force vector on a face is denoted by $f_e$. I now need to evaluate the integral $$\int_{\Gamma_e} N^T f_e \ d\Gamma_e $$ where $\Gamma_e$ represent any face. $N^T$ represent the shape function matrix given by $$\begin{bmatrix} \xi_1 & 0 & 0 & \cdots & \xi_4 & 0 & 0 \\ 0 & \xi_1 & 0 & \cdots & 0 & \xi_4 & 0 \\ 0 & 0 & \xi_1 & \cdots & 0 & 0 & \xi_4 \end{bmatrix} $$ To evaluate the surface integral one takes the trace of the matrix on the particular face. For example on face one, $\xi_1=0$. I now need to use a Gaussian Quadrature and not sure how to evaluate the transformed Jacobians. I am actually implementing the above procedure on a 10 node element. So actually confused on how to implement the surface quadratures. Another source of doubt is that in case of the linear element one can do a node lumping i,e distributing the surface load equally among all the nodes and that gives a very simple expression for the surface vector. Can this be also be done for a higher element?


2 Answers 2


The particular surface integral you want to calculate is basically a specific case of integrating a function over a surface defined in terms of two parametric coordinates. Lets first consider this general case. A fundamental relation from differential geometry is (using your notation)

$$ d{\bf \Gamma_e} = \left\{\begin{array}{c} \frac{\partial x}{\partial\xi_1} \\ \frac{\partial y}{\partial\xi_1} \\ \frac{\partial z}{\partial\xi_1} \end{array}\right\} \times \left\{\begin{array}{c} \frac{\partial x}{\partial\xi_2} \\ \frac{\partial y}{\partial\xi_2} \\ \frac{\partial z}{\partial\xi_2} \end{array}\right\} $$ where $\xi_1$ and $\xi_2$ are the two parametric coordinates of the surface. An integral of a function, $f$ over the surface would be written as

$$\int f d\Gamma_e = \int f(\xi_1, \xi_2) ||d{\bf \Gamma_e}||_2 d\xi_1 d\xi_2$$ where $||d{\bf \Gamma_e}||_2$ is just the magnitude of the vector.

Now lets consider the specific case of faces of the ten-node tetrahedron. In an isoparametric formulation, $x, y,$ and $z$ are written in terms of the shape functions, $N_I$, e.g. $$ x = N_I(\xi_1, \xi_2) x_I $$ where $ x_I$ are the values of $x$ at the nodes of the triangular face. The expressions for $N_I$ on a face can be obtained from the element shape functions by specializing them for a particular face, or equivalently, by using the shape functions for the six-node triangle. From these equations, you can easily evaluate the vector cross product above.

The most general way to actually perform the integration is by using the appropriate Gauss-Legendre rule for the triangle.

Finally, you asked about the difference between formally evaluating the consistent load integral and simply lumping equal values at the nodes. The equivalent nodal loads are quite different from the simple lumped values. For a simple, straight-sided element, the equivalent nodal loads are non-zero only at the mid-side nodes. If you have a fine enough mesh, the simple lumped approach will give you the correct answer but it is not optimum in the Galerkin-FEM sense.

  • $\begingroup$ Thanks Bill! I require a little more detailing. I have the shape functions of the 10 node tetrahedron. So if say nodes 123456 form a face do I just take the shape functions corresponding to these nodes (trace)? does $\Gamma_e$ appear like the normal vector of the face? Also if nodes 123 form the corner of the face I just take the shape functions corresponding to this say $\xi_1, \xi_2, \xi_3$ and use them as my parametric coordinates? Also how do I get the Jacobian of the transformation? Sorry I am confused by the numerous references I am using. Carlos Fellipa, Bathe, Reddy to name a few! $\endgroup$
    – kaush
    May 16, 2017 at 4:35
  • 2
    $\begingroup$ Yes, as I said in the post you can use the shape functions for the tet more or less as you describe. As you observed, one of the tet coordinates equals zero on a face. And on a face, one of the tri coordinates depends on the other two. Then you calculate the Jacobian from my first equation. Its "determinant" is just the magnitude of that vector. I don't recall seeing a discussion of this surface force calculation in any of the three references you mention; there is a brief discussion in Zienkiewicz. $\endgroup$ May 16, 2017 at 11:00
  • $\begingroup$ How do find the dependency of the tet coordinates? Also one more thing I observe is that if you take the Gaussian points as (0.5 0.5 0),(0 0.5 0.5) (0.5 0 0.5) and weight as 1/3 you get just the mid nodes contributing to the integral! Does that mean no contribution to the side nodes? Also my Jacobian is just the face normal is that right. When I integrate the unit function I get the surface area.. $\endgroup$
    – kaush
    Jun 9, 2017 at 7:53
  • $\begingroup$ The shape functions at the corner node $i$ are given as $\xi_i(2\xi_i-1)$. whereas the shape functions at the mid node of edge connecting $i$ and $j$ is given by $4 \xi_i \xi_j$. Hence at the Gauss points only one of the shape function in the definition of $N$ is non-zero. $\endgroup$
    – kaush
    Jun 9, 2017 at 8:14
  • $\begingroup$ @BillGreene Thank you so much! If it is not too much, could you kindly help provide which part of the Zienkiewicz book (The Finite Element Method, I assume) is discussing this? Vol 1 or 2? (I am using the 5th ed, but I think the structure largely remain unchanged) $\endgroup$
    – Shawn Wang
    Sep 27, 2018 at 17:08

You may please refer chapter two of Advanced Topics in Finite Element Analysis of Structures: With Mathematica and MATLAB Computations by M. Asghar Bhatti.

Section 2.4.5 of the book deals with the evaluation of surface integrals. A workout example of evaluating the surface pressure loads with Mathematica code is also given in the book. A Matlab code is also given in the book site for 8 node brick element.

  • 4
    $\begingroup$ This answer would benefit a lot from the information why this particular book/chapter are particularly useful to answer the OP's question, especially since it was not a simple reference-request question. $\endgroup$
    – Anton Menshov
    Jun 6, 2018 at 2:59

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.