[Vasilis]

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

[Jonas Holdeman]

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.

[Vasilis]

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

[Jonas Holdeman]

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

[Vasilis]

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!!

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

[Vasilis]

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?

[Ananda Himansu]

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.

[Vasilis]

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?

[Ananda Himansu]

No.

Well, of course, I do not infringe on your freedom to do anything you like, within the law. However, the value of your double integral would not come out correctly. First of all, your use of the symbol y' is ambiguous, given its use as a boundary-fitted coordinate in your previous post. I apologize for creating confusion of a similar sort by using the symbol f for the boundary curves when you had already used it for the integrand f(x,y) in your post (a fact which escaped my notice, as I was skimming your post in a hurry). Here, by y'(x) I assume you mean the derivative of y wrt x.

I see your new suggestion as an attempt to mimic the change of coordinates formula for a single integral, prompted perhaps by my inadvertent (and erroneous) conflation of the integrand with the boundary curve. This will not work, because it leaves the inner integral unperformed. However, your idea about change of variable is the same idea that underlies the mapping method, except that the latter must be done in 2D rather than 1D as you are suggesting. The y'(x) that you introduce corresponds to the Jacobian in 2D or higher dimensional change-of-variables in multiple integrals.

The correct method without (x',y') is as described in my previous post. Please reinterpret my previous post, retaining f(x,y) as the integrand, and this time using the notation y = gn(x) for the n'th boundary curve, so that the boundary curves are y = g1(x), y = g2(x), y = g3(x), etc. Basically, this uses x as the parametric coordinate for the boundary curves, which is not ideal as the gn can be multivalued and thus you lose a large part of the power of parametric representation of curves and domains. Also, you can see how awkward the integration becomes with variable limits.

Just as an aside, in the mapping method for a four-sided element, the boundary curves are parametrized by x' for two of the sides and by y' for the other two sides. In the mapping method, the boundary curves are transformed to "vector" or paired functions of parameters x' and y'. The domain of integration is given by the vector {x(x',y'), y(x',y')}, while the boundary curves take the form {x(-1,y'), y(-1,y')} for one boundary, {x(1,y'), y(1,y')} for another, {x(x',-1), y(x',-1)} for a third, and {x(x',1), y(x',1)} for the fourth. You can see how convenient this is for doing the double integral, and how elegant the mapping method is. The only price is to find the right mapping, and to compute the Jacobian.

To really understand why the mapping method is so attractive, try finding the area of an arbitrary quadrilateral (with straight sides, none of which is parallel to any other or to the x or y axes) by a double integral. The function f(x,y) is now simply unity, and the boundary curves g(x) are straight lines. Try it both with a mapping and without. If you are still confused, you should crack open a chapter on transformation of multiple integrals in a good calculus book.