When one solves numerically partial differential equations, using finite element methods, soon or later will encounter with the concept of isoparametric elements.
My goal here is not to explain the concept (I assume you have an idea), but to provide the solution to a problem that I faced up recently: I had to convert parametric or natural or local coordinates to global and vice versa.
A bit of background:
Every point in a domain is characterized by its coordinates (x, y ,z). These coordinates are called global. However in real applications it is not easy to work with the real or global coordinates, therefore it is common to use another coordinate system, where the location and the values of each point is a function of the location or values of the nodes of the element that contain the point. In other words we interpolate the value of a point from the values of the point corners.
For example, in a triangle where v1, v2, v3 are the values on the three corners, we can express the value of a point using the following formula:
vp=v1*N1(u,w)+v2*N2(u,w)+v3*N3(u,w)
where, N1, N2, N3 are called shape functions, u, w, are the parametric variables and vp is the interpolated value (it can be simply the coordinates or hydraulic conductivity, or even a vector e.g velocity)
The parametric variables usually span within [0 1] or [-1 1].
Below is a small matlab/octave code that may help you to understand the concept of parametric variables.
(In case you don't know, the shape functions of a 2D triangle are N1=u, N2=w, and N3=1-u-w)
hold on % create a window plot
[x y]=ginput(3); %pick manually the three corners of your triangle
z=[0.7; 0.1; 0.3];%Assign some random property on the nodes (eg. height)
p=[]; %create an empty variable
for u=0:0.1:1 %loop through the parametric variable
for w=0:0.1:1 %same hereif u+w<=1% for triangles u+w must be less than one
%otherwise the point is outside of the triangle
p=[p;u*x(1)+w*x(2)+(1-u-w)*x(3)... %find the global
u*y(1)+w*y(2)+(1-u-w)*y(3)... %coordinates for
u*z(1)+w*z(2)+(1-u-w)*z(3)]; %each u,w pair
end
end
end
plot3(x([1:3 1]), y([1:3 1]), [0 0 0 0]);
hold on;
plot3(x([1:3 1]), y([1:3 1]), z([1:3 1]),'r');
plot3(p(:,1), p(:,2), p(:,3),'o');
The above code will produce a red triangle in 3D, the projection of the triangle in the x-y plane and a series of points that lay on the plane
We see that we can interpolate the height value of every point within the triangle using very simple calculations.
So far we showed how you can find the global coordinates of a point if you know the local ones. Generally this is a straightforward task, yet the opposite can be cumbersome.
Bellow I will explain a general method where using the symbolic toolbox of matlab this becomes also a trivial task.
For a 2D triangle the problem is stated as follows, given the three corners of a triangle (x1,y1), (x2,y2), (x3,y3) and the coordinates of a point (xp,yp) find the parametric coordinates (u,w).
In case of triangles this is quite easy because the shape functions have relatively simple form. First we write the two equations:
In case of triangles this is quite easy because the shape functions have relatively simple form. First we write the two equations:
xp=x1*u+x2*w+x3*(1-u-w)
yp=y1*u+y2*w+y3*(1-u-w)
and then we solve a system of two equations with two unknowns (u and w). This is actually something that matlab symbolic toolbox is very good at.
Lets first create few symbolic variables
>> syms x1 x2 x3 y1 y2 y3 xp yp u w
>> eq1 = x1*u+x2*w+x3*(1-u-w) - xp
>> eq2 = y1*u+y2*w+y3*(1-u-w) - yp
Finally, all we have to do is to tell matlab to solve this system for us
>> S=solve(eq1==0,eq2==0,u,w);
The outcome of the last command is a structured variable with fields u and w, which contains the solutions for u and w.
S.u=(x2*y3 - x3*y2 - x2*yp + xp*y2 + x3*yp - xp*y3)/(x1*y2 - x2*y1 - x1*y3 + x3*y1 + x2*y3 - x3*y2)
S.w=-(x1*y3 - x3*y1 - x1*yp + xp*y1 + x3*yp - xp*y3)/(x1*y2 - x2*y1 - x1*y3 + x3*y1 + x2*y3 - x3*y2)
Of course for the 2D triangle one can derive these formulas using a different approach.
