Shape function gradients with 3D boudary elements
 

Hello all,

I am developing a 3D boundary element program where I have only surface elements (4-node quads and 3-node tris), but that are oriented in 3D space. Likewise, I want to compute x,y,z gradients on each of these faces.

I have a nodal solution (phi), and I want to compute dphi/dx, dphi/dy, dphi/dz.

I have successfully done this for quads, but cannot seem to figure it out for triangles.

For the quad-4 with linear shape functions:
N = [1/4*(1+xi)*(1+eta) 1/4(1+xi)*(1-eta) 1/4(1-xi)*(1-eta) 1/4(1-xi)*(1-eta)]
I can compute the differential operators in x,y,z (dNdx, dNdy, and dNdz) by simply computing 2 tangential vectors (dxdxi, dxdeta, dydxi, dydeta, dzdxi, dzdeta):
e1 = dNdxi*[ex,ey,ez]
and
e2 = dNdeta*[ex,ey,ez]

Taking the cross product of these gives me my surface gradient:
e3 = cross(e1,e2)

Then using these 3 vectors as my jacobian:
Jac = [e1; e2; e3];

I can invert my jacobian to compute dxidx,dxidy,dxidz, detadx,detady,detadz, and finally post multiply with dNdxi and dNdeta again to compute dNdx, dNdy, dNdz:

dNdx = Jac(1,1)^-1*dNdxi + Jac(1,2)^-1*dNdeta
dNdy = Jac(2,1)^-1*dNdxi + Jac(2,2)^-1*dNdeta
dNdz = Jac(3,1)^-1*dNdxi + Jac(3,2)^-1*dNdeta

These work well, and I get the correct gradients.

But now for the 3-node triangle, it's not so straight forward. I have formulated the TRI-3 with 3 local coordinates (xi, eta, and zeta), where my shape functions are just:
N = [xi eta zeta],
and so:
dN = I (identity).

How can I compute dNdx, dNdy, and dNdz for the linear 3-node triangle?

Any help is appreciated, and thank you in advance.

I am not sure about my answer, so sorry if I have not understood your question.

Consider a 3D tethraedron (4-nodes) where you define linear shape function Ni(x,y,z) (i=1,..4).

The Grad operator is: i d/dx + j d/dy + k d/dz (in bold the unit vectors in the Cartesian reference system).
Now what is the problem in computing analytically the gradient?

The quad element has curvature; the tri element with linear shape functions does _not_ have curvature! I think this should explain your problem.

adrin

Hi, thanks for the reply.

I am referring to a 3D boundary element, not a 3D volume element. So image that I want to compute a gradient in 3 dimensions, on the face the tetrahedron.

This is what I am trying to do.

Here, you cannot use standard 2D formulations, but since it is a boundary element, standard 3D volume formulations also do not apply.

Hi, thanks for the reply.

I agree that a quad could have a torsional curvature to it, and the triangle cannot represent this, however, this has no bearing on computing the derivatives over the surface. I agree that the normal would vary across the panel, and this could be computed using the formulation in my original post, where xi and eta would be evaluated at various points (say guass points), and the normal would be calculated to be different across the element.

However, leaving this issue out of it (say, for example, that I do not allow any torsion in any quads), I still cannot derive how to compute a discrete exterior gradient over the triangular element.

I am confused in terms of the problem you want to solve...
think a simple example of boundary element method such as the panel method on airfoil. You MUST fix the normal component of the velocity that directly imply a condition on the normal derivative. So I think that what you want to compute is actually a BC dependent issue.

However, it the shape function are linear, any derivative produce a constant value

I'm having a bit of a hard time following your formulations/notations (even for the quad element). N is a 4-element vector, so how does dN/dxi*[ex,ey,ez] lead to e1, which appears to be a 3-element vector?

For the Tri element, zeta is not independent; it should be zeta = 1 - xi - eta. Also, I don't get how you arrived at dN = I (part of it has to do with dN vs dN/d(something) and the other part is with I itself; is it a vector, tensor, ...?)

adrin

Just take a triangle arbitrarily oriented in space. Each corner of the triangle has a node, and with that node, x,y,z coordinates, and some value phi. eg. Node (i) has the coordinates x_i,y_i,z_i, and phi value phi_i.

I want to linearly interpolate over the triangular element, so I use the interpolation functions N = [xi eta zeta].

where xi, eta, and zeta are the local coordinates aligned with each corner->side direction.

Taking the derivative of N with respect to xi, eta, and zeta, you get I = [1 0 0; 0 1 0; 0 0 1].

I now want to linearly interpolate the gradient of phi across the triangular element: dphi/dX. where X = [x,y,z]. eg. dphi/dX = [dphi/dx dphi/dy dphi/dz].

This is what I am trying to compute.

Thanks for your interest! Sorry, let me clarify.

For the QUAD4, N is [0.25*(1-xi)*(1+eta) 0.25*(1+xi)*(1+eta) 0.25*(1+xi)*(1-eta) 0.25*(1-xi)*(1-eta)]

Taking the derivative with respect to xi and eta, dN/dxi = [-0.25*(1+eta) 0.25*(1+eta) 0.25*(1-eta) -0.25*(1-eta)] and dN/deta = [0.25*(1-xi) 0.25*(1+xi) -0.25*(1+xi) -0.25*(1-xi)]

letting dN = [dN/dxi ; dN/deta], dN will be a 2 x 4 matrix.

e1 = dN/dxi*[ex ey ez] is a 1x3 row vector because ex, ey, and ez are 4 component column vectors, and dNdxi is a 1x4 row vector.

Likewise, e2 = dN/deta*[ex ey ez], is a 1x3 row vector.

e3 = e1 x e2, is a 1x3 row vector.

J = [e1;e2;e3], is the 3x3 Jacobian matrix.
Then by inverting the Jacobian, the gradients of the shape functions are calculated:

dNdx = J^-1(1,1:2)*dN
dNdy = J^-1(2,1:2)*dN
dNdz = J^-1(3,1:2)*dN


For the 3 node triangle, see my above response to @FMDenaro how I got dN = I.I am now looking into using 2 independent local coordinates instead of 3 to see if I can analogously compute the gradient like with the quad.

I am not sure about my proposal, however check it.
in the (x,y,z) the Grad phi = i dphi/dx + j dphi/dy + k dphi/dz that can be equate to the expression of the gradient in the local (orthogonal) reference system on the triangle. That means for example that from the equality of the expression you can compute dphi/dx by scalar product with i and so on for the other derivatives. In conclusion you need only the angles of orientation

Note, formulating the triangle element as 2 local coordinates, rather than three, you can just use the same theory as the quad above.

Set zeta = 1-xi-eta and proceed.

using normalized areal coordinates that is correct