2D Quadrilateral Interpolation: Barycentric Coordinates

Interpolation on a general quadrilateral

When dealing with general quadrilaterals, the explicit expression of the shape functions is not easy to found. But we have the possibility of computing the value of such functions on a point p using the barycentric coordinates for quadrilaterals.

Similar to the triangular interpolation, now consider a point

, in a quadrilateral defined by four vertices

numbered counterclokwise as usual(see figure). One can define the barycentric coordinades of this point which are equivalent to the value of the 2D shape functions

defined for this quadrilateral

satisfying

This way, one can express any point using four coordinates (barycentric)

. In particular

The barycentric coordinates satisfy

. Moreover all

unless the point is outside of the quadrilater.

Thus, the coordinates of the point

To compute the values of

we re-write the previous equation as

where

Therefore,

Now to compute these two parameters, we rename the terms in the previous equation as:

where

Finally we apply Newton's method to obtain

Compute Barycentric Coordinates on a general quadrilateral

Given a quadrilateral defined by vertices

, compute the barycentric coordinates of point

Verify that

% Data values

p=[3 2];

v1=[0 0];

v2=[5 -1];

v3=[4 5];

v4=[1 4];

% define the previous notation

a=(v1-p);

b=v2-v1;

c=v4-v1;

d=v1-v2-v4+v3;

% we use (x,y) instead of (lambda, mu) for easy code

% solution is supposed to have 0 <= x <=1 and 0 <= y <=1

x=0.5; y=0.5; %initial point

dx=0.1*ones(2,1);

iter=0;

tol=1.e-12;

while(norm(dx) > tol && iter < 20 )

%apply Newton-Rapson method to solve f(x,y)=0

f = a + b*x + c*y + d*x*y;

% Newton: x_{n+1} = x_n - (Df^-1)*f

% or equivalently denoting dx = x_{n+1}-x_n

% Newton: Df*dx=-f

Df(:,1)=(b + d*y); %df/dx

Df(:,2)=(c + d*x); %df/dy

bb=-f'; %independent term

dx=Df\bb;

x=x+dx(1);

y=y+dx(2);

iter=iter+1;

if (norm([x,y]) > 10), iter=20; end %non convergent: just to save time

end

if (iter < 20)

lambda=x;

mu=y;

alpha1=(1-mu)*(1-lambda);

alpha2=lambda*(1-mu);

alpha3=mu*lambda;

alpha4=(1-lambda)*mu;

alphas=[alpha1,alpha2,alpha3,alpha4]

else

alphas=[-1,-1,-1,-1]; %wrong values

end

% obtained point

q=alpha1*v1+alpha2*v2+alpha3*v3+alpha4*v4;

p-q %we must recover the same point

% Plot the results

plotRectangle(v1,v2,v3,v4);

hold on;

plot(p(:,1),p(:,2),'o');

plot(q(:,1),q(:,2),'*');

hold off;

Exercise 1:

From the previous code implement the function baryCoordQuad (the Quadrilateral version of the baryCoord function).

function [alphas,isInside]=baryCoordQuad(vertices,point)

Exercise 2:

From the mesh file meshTwoHolesQuad.m, find the quadrilateral where the point p = [39, 7] belongs to.

Sol:

i = 404

vertices =

41.9240 8.7481

37.5060 8.6827

36.5394 4.7543

41.0196 4.8224

Hint:

Be carefull with the final results, you must check that isInside = 1.

In the case you ask for a point that is not inside the meshed domain, then, no element meet the condition of isInside = 0. You can cheek this, using point p = [0,0] for the previous example (is the center of one of the holes)

Exercise 3:

Using the same mesh and point as in the previous exercises, consider the nodal temperature defined by the function:

Compute the interpolated temperature at p = [39, 7]:

Sol:

clear all
% plot a mesh
eval('meshTwoHolesQuad'); %load data from the file mesHole.m: nodes and elements
numElem = size(elem,1);  %number of mesh elements
numNodes = size(nodes,1); %number of mesh nodes
p=[39,7];
% evaluate temperature function f(x,y) = 100/(x^2+y^2);
x=nodes(:,1);
y=nodes(:,2); 
uTemp=100./(x.^2+y.^2);

 
%% Version 1
%
% Look for the element the point p belongs to:
%
for e=1:numElem % for each element
    n1 = elem(e,1); % number of the 1st node of the element 
    n2 = elem(e,2); % number of the 2nd node of the element 
    n3 = elem(e,3); % number of the 3th node of the element
    n4 = elem(e,3); % number of the 4th node of the element
    v1=nodes(n1,:); % coordinates of the 1st node of the element 
    v2=nodes(n2,:); % coordinates of the 2nd node of the element 
    v3=nodes(n3,:); % coordinates of the 3rd node of the element 
    v4=nodes(n4,:); % coordinates of the 4th node of the element 
    vertices=[v1; v2; v3; v4];
    [alphas, isInside]=baryCoordQuad(vertices, p);
    if (isInside == 1) % Check if p is inside, compute the results and break the FOR loop
        elemP = e;             %save the present element
        nodElemP = [n1,n2,n3,n4]; %save the number nodes
        alphasP = alphas;      %save the barycentric coordinates
        break   %exit from the for loop
    end 
end
% use barycentric coordinates to interpolate vertices values
uPointInterp = alphasP(1)*uTemp(n1) + alphasP(2)*uTemp(n2) + alphasP(3)*uTemp(n3)+alphasP(4)*uTemp(n4);
fprintf('The point belongs to element %d;\n', elemP);
fprintf('Their node numbers are %d, %d, and %d.\n', nodElemP);
fprintf('The interpolated value is uPInterpol = %6.4f\n',uPointInterp);
 

%% Version 2: using vector notation
%
% Look for the element the point p belongs to:
%
for e = 1:numElem     % for each element
    vertices = nodes(elem(e,:),:); % vertices of element e
    if size(elem,2)==3 %triangles
        [alphas, isInside] = baryCoord(vertices, p); %Barycentric coords
    elseif size(elem,2)==4 %Quadrilaterals
        [alphas, isInside] = baryCoordQuad(vertices, p); %Barycentric coords
    end
    if (isInside == 1) % Check if p is inside, compute the results and break the FOR loop
        elemP = e;             %save the present element
        nodElemP= elem(e,:); %save the number node
        alphasP = alphas;      %save the barycentric coordinates
        break   %exit from the for loop
    end 
end
% use barycentric coordinates to interpolate vertices values
uPInterpol = dot(alphasP,uTemp(nodElemP))
fprintf('The point belongs to element %d;\n', elemP);
fprintf('Their node numbers are %d, %d, and %d.\n', nodElemP);
fprintf('The interpolated value is uPInterpol = %6.4f\n',uPointInterp);

(c)Numerical Factory 2020