# Need help to use Hermite basis functions

 Register Blogs Members List Search Today's Posts Mark Forums Read

 October 22, 2009, 10:20 Need help to use Hermite basis functions #1 New Member   Join Date: Apr 2009 Posts: 17 Rep Power: 8 Dear all, I am trying to solve a fairly simple test case problem with Hermite basis function and FEM. Lets suppose that F_i is the set of basis functions, and this set is comprised of two sets: H_i and S_i. H_i are the basis functions whose value is equal to one at the nodes, and S_i are basis functions whose derivative is equal to one at the nodes. If u_e is the solution, then it is approximated as u_e=Sum u_i * F_I I will calculate the integral with the Gauss quadrature method, thus I will use an iso-mapping to go from the parent to the computational element. x=Sum x_i * H_i + Sum xx_i * S_i, where xx_i is the derivative of x_i w.r.t. ksi (ksi takes value in the range -1,1) I have read on the Internet that in order to compute xx_i, I need to consider a perturbation in the neighborhood of x_i, and then calculate the xx_i. I decided to do something much simpler. I considered a new set of basis functions f_i (so that their value is equal to one at the nodes), and I used this set to compute xx_i=Sum x_i * fx_i, where fx_i is the derivative of f_i w.r.t ksi. Do all this make sense? Of course the code does not converge, and I am trying to understand what the problem might be. I do not think that it is in the way I am calculating xx_i ... but, could it be? Thank you for your time.. EDIT: There are only two nodes per element

 October 23, 2009, 10:27 #2 Member   Jonas T. Holdeman, Jr. Join Date: Mar 2009 Location: Knoxville, Tennessee Posts: 90 Rep Power: 8 Perhaps I don't understand what you are saying, but it seems to me that you are not really using Hermite functions/elements correctly. Let me reuse your notation. Let F_i be the set of Hermite functions. These functions have two parts, H_i which take the value 1 at node i and 0 at all other nodes with derivatives taking values 0 at all nodes, and functions S_i which vanish at all nodes, with first derivative vanishing at all nodes except node i where the derivative is 1. Then you would write u_e = Sum(u_i H_i+ux_i S_i) If you write these as vectors, i.e. let U_i=[u_i, ux_i]^T and F_i=[H_i, S_i]^T (where ^T indicates the transpose), then we could write the vector equation, u_e = Sum(F_i^T U_i) = Sum(U_i^T F_i) where we are summing over dot products which gives a scalar for u_e. You could use the Hermite functions for an isoparametric mapping in which case you would specify x_i and xx_i at each node (or the vector X_i = [x_i, xx_i]^T), in which case the xx_i are specified as part of the mesh generation. It is much simpler though to use a linear parametric mapping. This would be equivalent to using elements with straight sides in 2D. I have used modified Hermite functions/elements quite successfully for incompressible fluids in 2D, where the h_i represent the stream function and the degrees of freedom of s_i are components of the velocity (I am changing notation here). Then if we write S_i^T = curl([h_i, s_i]) we have, U_e = Sum(S_i^T U_i) where the components of U_i are the stream function and velocity components at node i.

 October 23, 2009, 10:52 #3 New Member   Join Date: Apr 2009 Posts: 17 Rep Power: 8 Well, I need to clarify one thing. Based on you notation, the u_i and ux_i are both unknowns that we need to solve for. This means that if each element has 2 nodes, we have 4 unknowns (2 for u_i and 2 for ux_i). Now, the iso-mapping has the form: x=Sum (x_i * H_i + xx_i * S_i), and xx_i is the equivalent to ux_i, but it is not an unknown, and it has to be computed, based on the value of x_i. You wrote that xx_i are specified as part of the mesh generation, can you give me an example? If the mesh is a straight line, are the xx_i=0? Things are pretty complicated if I have an oblique mesh, however!! I managed to solve the problem (both solution and slope are correct) by not using a parent element, and thus I did not need to use any kind of mapping. Although, this method is straight forward, I would like to be able to solve the problem with the use of a parent element. A parent element is very convenient in moving boundaries problems. Now, in the past couple of hours, I 've being trying different schemes for x. If I omit the term xx_i * S_i (i.e. xx_i=0, since the mesh is a straight line), I do solve the problem and I do get the correct solution, but the slope is totally wrong. Similar results I get if I omit the x_i * H_i term. Could you elaborate on "linear parametric mapping" and elements with straight sides? I intend to use 3 nodes per element (1D) or 9 nodes per element (2D). Thank you

 October 23, 2009, 13:34 #4 Member   Jonas T. Holdeman, Jr. Join Date: Mar 2009 Location: Knoxville, Tennessee Posts: 90 Rep Power: 8 You wrote that xx_i are specified as part of the mesh generation, can you give me an example? I assume that you are using a graded mesh. If your nodal coordinates are given by some function f so that x_i = f(i), then xx_i can be found from the slope of f at node i. If only the x_i are specified, then plot x_i versus i. If these points are connected by straight line segments, then this would correspond to a linear mapping. If you insist on (cubic) Hermite mapping, then xx_i^e on each interval would be the slope of the straight line segment. Note that you cannot assign xx_i to each node because the slope is discontinuous there. However, if you interpolate the points (or groups of points) with a smooth function/curve, the xx_i can be found from the slope of the interpolation. This is all done without reference to the dependent problem variables, which is why I referred to it as part of mesh generation. If the mesh is a straight line, are the xx_i=0? No, the xx_i would be the slope of the straight line. If you set xx_0=0 at the nodes, the curve of the mapping would be a wiggly stairstep function. Could you elaborate on "linear parametric mapping" and elements with straight sides? In 2D I should have said bilinear mapping of Hermite elements. If the local variables of the parent element N_i are (q,r) on [-1,1]x[-1,1], then we would have x(q,r) = Sum(x_i*N_i(q,r)), and similarly for y. On element boundaries we have q=+/- 1 or r=+/- 1, and x and y are linear in the remaining variable and these represent straight lines in the global coordinate plane. Using simple cubic or bicubic mapping would, in general, give curved sides. I intend to use 3 nodes per element (1D) or 9 nodes per element (2D). Using u_i and ux_i at each of nine nodes would imply quintic elements. I have done this, but I prefer to use four-node quintics in 2D. If you want to see some results of incompressible fluid flow using Hermite elements with bilinear mapping, look at the paper at http://webpages.charter.net/jtholdeman/RevIJNMFart1.pdf

October 26, 2009, 07:12
#5
New Member

Join Date: Apr 2009
Posts: 17
Rep Power: 8
Quote:
 Originally Posted by Jonas Holdeman In 2D I should have said bilinear mapping of Hermite elements. If the local variables of the parent element N_i are (q,r) on [-1,1]x[-1,1], then we would have x(q,r) = Sum(x_i*N_i(q,r)), and similarly for y. On element boundaries we have q=+/- 1 or r=+/- 1, and x and y are linear in the remaining variable and these represent straight lines in the global coordinate plane. Using simple cubic or bicubic mapping would, in general, give curved sides.
This type of mapping is similar to the one I wrote. The only difference is that you omitted the term xx_i * S_i. I tried to use this mapping, and although I got the correct solution, the slope is totally wrong!!

Quote:
 If the mesh is a straight line, are the xx_i=0? No, the xx_i would be the slope of the straight line. If you set xx_0=0 at the nodes, the curve of the mapping would be a wiggly stairstep function.
Yes, you are correct, but since I am talking about 1D problem, the mesh is a straight line which is parallel to the x-axis. So, xx_i=0.

Although I hate saying this, I will stick to the global basis functions. Thus, I do not have to use any mapping, and until now it is working.

Thank you Jonas for all your help

 October 29, 2009, 09:08 #6 New Member   Join Date: Apr 2009 Posts: 17 Rep Power: 8 Here is an interesting question: When one uses any kind of mapping, one can transform the double integral from dxdy to |J|dx'dy', where |J| is the determinant of the Jacobean matrix, dx' and dy' are the independent variables in the parent element. The lower limits of the integrals are -1 and the upper limits are 1, i.e., int_a_b { int_c_d} f(x,y)dx }dy ==> int_-1_1 {int_-1_1 f(x',y')dx' }dy' Now, if I do not use a mapping, how would I calculate a double integral? If I know the shape of the element, I could find the expression that relates y to x, and convert the double integral to a single integral. But, if I do not know a priori the shape of the element, how do I compute the double integral dxdy?

 October 29, 2009, 09:32 #7 New Member   Ananda Himansu Join Date: Apr 2009 Location: Cleveland, Ohio, USA Posts: 17 Rep Power: 8 Mapping is exactly what turns the boundary curves in the (x, y) plane to axes-parallel straight line segments in the (x', y') plane. If you do not wish to map the element to a rectangle, then you have to do the double integral as an iterated (nested) integral, with variable limits of integration on the inner integral. Let us say you know the boundary curves of the element in the form y = f(x). Then, by inverting these boundary curves and minimizing and maximizing, you have to find the limits of integration xmin and xmax. The limits on the outer integral will then by xmin and xmax. The inner integral will have variable limits ymin and ymax which are a function of x, and will depend on the boundary curves f1(x), f2(x), f3(x), ... Note that this can approach can cause problems if the tangent to a boundary curve is vertical at any point or segment of it.

October 29, 2009, 09:56
#8
New Member

Join Date: Apr 2009
Posts: 17
Rep Power: 8
Quote:
 Originally Posted by Ananda Himansu Let us say you know the boundary curves of the element in the form y = f(x).
If I know this boundary curve, cann't I substitute y by f(x) in the integrand, multiply the integrand by y'(x) and then integrate over x?