interpolation_on_mesh.m
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% This file is part of NS2DDV. %
% %
% Copyright(C) 2011-2018 C. Calgaro (caterina.calgaro@math.univ-lille1.fr) %
% E. Creusé (emmanuel.creuse@math.univ-lille1.fr) %
% T. Goudon (thierry.goudon@inria.fr) %
% A. Mouton (alexandre.mouton@math.univ-lille1.fr) %
% %
% NS2DDV is free software: you can redistribute it and/or modify it under the terms %
% of the GNU General Public License as published by the Free Software Foundation, %
% either version 3 of the License, or (at your option) any later version. %
% %
% NS2DDV is distributed in the hope that it will be useful, but WITHOUT ANY %
% WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A %
% PARTICULAR PURPOSE. See the GNU General Public License for more details. %
% %
% You should have received a copy of the GNU General Public License along with %
% NS2DDV. If not, see <http://www.gnu.org/licenses/>. %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function[unew] = interpolation_on_mesh(uold, METHOD, MESHOLD, MESHNEW, PARAMETERS)
% INPUT
% uold Data discretized on nodes of old mesh
% METHOD FE method associated to the discretization of uold
% MESHOLD Old mesh
% MESHNEW New mesh
% PARAMETERS List of simulation parameters
%
% OUTPUT
% unew Data interpolated on nodes of new mesh
switch METHOD
case {'P1'}
pold = MESHOLD.p1;
npold = MESHOLD.np1;
told = MESHOLD.t1;
ntold = MESHOLD.nt1;
pnew = MESHNEW.p1;
npnew = MESHNEW.np1;
tnew = MESHNEW.t1;
ntnew = MESHNEW.nt1;
nDOF = 3;
case {'P2'}
pold = MESHOLD.p2;
npold = MESHOLD.np2;
told = MESHOLD.t2;
ntold = MESHOLD.nt2;
pnew = MESHNEW.p2;
npnew = MESHNEW.np2;
tnew = MESHNEW.t2;
ntnew = MESHNEW.nt2;
nDOF = 6;
case {'P1B'}
pold = MESHOLD.p1b;
npold = MESHOLD.np1b;
told = MESHOLD.t1b;
ntold = MESHOLD.nt1b;
pnew = MESHNEW.p1b;
npnew = MESHNEW.np1b;
tnew = MESHNEW.t1b;
ntnew = MESHNEW.nt1b;
nDOF = 4;
otherwise
error('Unknown FE method');
end
cprec = 1e-10;
loctab = zeros(npnew,1);
unew = zeros(npnew,size(uold,2));
switch PARAMETERS.PARALLELIZATION.TOOLBOX
case {'PCT'}
parfor ipnew=1:npnew
% We look for the triangle in the old mesh where the new node is located
P = pnew(:,ipnew);
for itold=1:ntold
% We compute the barycentric coordinates of the new node in the old triangle
p1 = pold(:,told(1,itold));
p2 = pold(:,told(2,itold));
p3 = pold(:,told(3,itold));
lambda1 = ((p2(2)-p3(2))*(P(1)-p3(1)) + (p3(1)-p2(1))*(P(2)-p3(2))) / ...
((p2(2)-p3(2))*(p1(1)-p3(1)) + (p3(1)-p2(1))*(p1(2)-p3(2)));
lambda2 = ((p3(2)-p1(2))*(P(1)-p3(1)) + (p1(1)-p3(1))*(P(2)-p3(2))) / ...
((p2(2)-p3(2))*(p1(1)-p3(1)) + (p3(1)-p2(1))*(p1(2)-p3(2)));
lambda3 = 1. - lambda1 - lambda2;
in = ((max(max([lambda1,lambda2,lambda3])) <= 1+cprec) && ...
(min(min([lambda1,lambda2,lambda3])) >= -cprec));
if in
% The current old triangle will be used for interpolating on the new node
loctab(ipnew) = itold;
break;
end
end
if (loctab(ipnew) == 0)
error(sprintf('Node %d has not been localized (%g,%g)', ipnew, P(1), P(2)));
else
% At this point, we have found an old triangle that will be used to interpolate on the current
% new node
% We evaluate the old DOF functions on the new node
% func_base doit être évalué dans le triangle de référence (utiliser phi_T^-1)
itoldloc = loctab(ipnew);
p1 = pold(:,told(1,itoldloc));
p2 = pold(:,told(2,itoldloc));
p3 = pold(:,told(3,itoldloc));
M = [p3(2)-p1(2), p1(1)-p3(1); p1(2)-p2(2), p2(1)-p1(1)];
d = (p2(1)-p1(1))*(p3(2)-p1(2)) - (p2(2)-p1(2))*(p3(1)-p1(1));
Pref = M*(P-p1)/d;
oldphiloc = func_base(METHOD, 0, Pref');
uoldloc = uold(told(1:nDOF,itoldloc),:);
unew(ipnew,:) = oldphiloc*uoldloc;
end
end
case {'SERIAL'}
for ipnew=1:npnew
% We look for the triangle in the old mesh where the new node is located
P = pnew(:,ipnew);
for itold=1:ntold
% We compute the barycentric coordinates of the new node in the old triangle
p1 = pold(:,told(1,itold));
p2 = pold(:,told(2,itold));
p3 = pold(:,told(3,itold));
lambda1 = ((p2(2)-p3(2))*(P(1)-p3(1)) + (p3(1)-p2(1))*(P(2)-p3(2))) / ...
((p2(2)-p3(2))*(p1(1)-p3(1)) + (p3(1)-p2(1))*(p1(2)-p3(2)));
lambda2 = ((p3(2)-p1(2))*(P(1)-p3(1)) + (p1(1)-p3(1))*(P(2)-p3(2))) / ...
((p2(2)-p3(2))*(p1(1)-p3(1)) + (p3(1)-p2(1))*(p1(2)-p3(2)));
lambda3 = 1. - lambda1 - lambda2;
in = ((max(max([lambda1,lambda2,lambda3])) <= 1+cprec) && ...
(min(min([lambda1,lambda2,lambda3])) >= -cprec));
if in
% The current old triangle will be used for interpolating on the new node
loctab(ipnew) = itold;
break;
end
end
if (loctab(ipnew) == 0)
error(sprintf('Node %d has not been localized (%g,%g)', ipnew, P(1), P(2)));
else
% At this point, we have found an old triangle that will be used to interpolate on the current
% new node
% We evaluate the old DOF functions on the new node
% func_base doit être évalué dans le triangle de référence (utiliser phi_T^-1)
itoldloc = loctab(ipnew);
p1 = pold(:,told(1,itoldloc));
p2 = pold(:,told(2,itoldloc));
p3 = pold(:,told(3,itoldloc));
M = [p3(2)-p1(2), p1(1)-p3(1); p1(2)-p2(2), p2(1)-p1(1)];
d = (p2(1)-p1(1))*(p3(2)-p1(2)) - (p2(2)-p1(2))*(p3(1)-p1(1));
Pref = M*(P-p1)/d;
oldphiloc = func_base(METHOD, 0, Pref');
uoldloc = uold(told(1:nDOF,itoldloc),:);
unew(ipnew,:) = oldphiloc*uoldloc;
end
end
otherwise
error('Unknown parallelization technique');
end
end
