Thursday, August 23, 2012

Inverse transformation of parametric elements

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 here
      if 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:
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.


Now lets try to do the same for the quadrilateral element:
Quadrilateral elements are any kind of  four point convex polygons.
The parametric variables u, w in that case span from -1 to 1.
Similar to triangle,  the value of any point is equal to:
vp=N1*v1+N2*v2+N3*v3+N4*v4
The shape functions for the quadrilateral are:
N1=0.25.*(1-u).*(1-w)
N2=0.25.*(1+u).*(1-w)
N3=0.25.*(1+u).*(1+w)
N4=0.25.*(1-u).*(1+w)
Lets first create a plot similar to the one above for the triangle
hold on % create a window plot
[x y]=ginput(4); %pick manually the four corners of the quad
z=[0.7; 0.1; 0.3; 0.5];%Assign some random property on the nodes
p=[]; %create an empty variable
N1=inline('0.25*(1-u)*(1-w)'); create the shape functions
N2=inline('0.25*(1+u)*(1-w)');
N3=inline('0.25*(1+u)*(1+w)');
N4=inline('0.25*(1-u)*(1+w)');  
for u=-1:0.1:1 %loop through the parametric variable
  for w=-1:0.1:1 %same here
     p=[p;N1(u,w)*x(1)+N2(u,w)*x(2)+N3(u,w)*x(3)+N4(u,w)*x(4)...            N1(u,w)*y(1)+N2(u,w)*y(2)+N3(u,w)*y(3)+N4(u,w)*y(4)... 
          N1(u,w)*z(1)+N2(u,w)*z(2)+N3(u,w)*z(3)+N4(u,w)*z(4)];
  end
end
plot3(x([1:4 1]), y([1:4 1]), [0 0 0 0 0]);
hold on;
plot3(x([1:4 1]), y([1:4 1]), z([1:4 1]),'r'); 
plot3(p(:,1), p(:,2), p(:,3),'o');

You can see that for the quadrilateral the interpolated points do not lay on the same plane.











To express the coordinates of a point as a function of the coordinates of the corner points we should write the following two equations:
xp=0.25.*(1-u).*(1-w)*x1+0.25.*(1+u).*(1-w)*x2+0.25.*(1+u).*(1+w)*x3+0.25.*(1-u).*(1+w)*x4
yp=0.25.*(1-u).*(1-w)*y1+0.25.*(1+u).*(1-w)*y2+0.25.*(1+u).*(1+w)*y3+0.25.*(1-u).*(1+w)*y4
If we try to solve this system we will end up with huge expressions. Personally I find it impossible for me to carry out by hand all the calculations without making mistakes.
So I always try to do such things with matlab
First lets define few symbolic variables
syms x1 x2 x3 x4 y1 y2 y3 y4 xp yp u w
N1=0.25.*(1-u).*(1-w); %define shape functions
N2=0.25.*(1+u).*(1-w);%since u is symbolic variable N will be symbolic expressions
N3=0.25.*(1+u).*(1+w);% and not inline functions
N4=0.25.*(1-u).*(1+w);
eq1=x1*N1+x2*N2+x3*N3+x4*N4-xp;
eq2=y1*N1+y2*N2+y3*N3+y4*N4-yp;
All we have to do is to tell matlab to try to solve the system.
S=solve(eq1==0,eq2==0,u,w);
You will notice that there are two solutions. and you have to try which ones gives the desired outcome( in my case it was u=S.u(1) and w=S.w(1)). You may also notice that both expressions of u and w are extremely lengthy!

In the code below, first you pick four corners (make sure the polygon is convex).
Next you pick two points in the global elements. The code plots few points along the line defined by the to points and the transforms the global coordinates to local using the symbolic solutions. (Note that the evaluation of the symbolic expression can be significantly slower, therefore in real application it is highly recommended to copy paste the expression in the code).  Finally the few points along the line are plotted in the local system.

[x y]=ginput(4);
x=100*x;
y=100*y;
subplot(1,2,1);plot(x([1:4 1]),y([1:4 1]))
hold on
subplot(1,2,2);plot([-1 -1 1 1 -1],[-1 1 1 -1 -1]);
hold on
subplot(1,2,1);title('Hit enter to exit')
while 1
    p1=[];p2=[];
    p1=ginput(1);
   
    if isempty(p1)
        break
    end
    p2=ginput(2);
    p(:,1)=p1(1,1)+(p2(1,1)-p1(1,1)).*[0:0.1:1]';
    p(:,2)=p1(1,2)+(p2(1,2)-p1(1,2)).*[0:0.1:1]';
    subplot(1,2,1);plot(p(:,1),p(:,2),'*')
    for i=1:length(p)
        u(i,1)=subs(S.u(1),[x1 x2 x3 x4 y1 y2 y3 y4 xp yp],[x' y' p(i,1) p(i,2)]);
        w(i,1)=subs(S.w(1),[x1 x2 x3 x4 y1 y2 y3 y4 xp yp],[x' y' p(i,1) p(i,2)]);
    end
    subplot(1,2,2);plot(u,w,'*')
end
            In the example the letter A looks upside-down in the local space because the first point of the polygon is the (46, 92) top corner. The first corner corresponds always to the corner with (-1, -1) in the local system.
Global coordinates                         Natural coordinates
                










             


No comments: