pmicro2tmacro.m

function pmic2tmac = pmicro2tmacro(glob,model)
%function pmic2tmac = pmicro2tmacro(glob,model)
%
% function defining a vector pmic2tmac containing the information
% which point of the vector glob lies in which triangle of the
% macrogrid.
% pmic2tmac(5) = 7 means that the 5th point in glob (glob(5,:) 
% lies in macro-triangle nr 7
%
% input:
% grid : triagrid
% glob : n-by-2 vector with the coordinates of the points
%
% output:
% mic2mac : n-by-1 vector 
%
% see also function micro2macro_map, which does the same fot microtriangles
% instead of single points
%
% Oliver Zeeb, 08.02.11

grid=triagrid(model.pmacro',model.tmacro',[]);


pmic2tmac = zeros(size(glob,1),1);

pmacro_x = grid.X;
pmacro_y = grid.Y;
tmacro = grid.VI;
nr_macro_tri =grid.nelements; %nr of macro-triangles

pmic_x = glob(:,1);
pmic_y = glob(:,2);

%dummy matrices for the affine transfomation from original macro triangle
%to reference triangle
C=zeros(2,nr_macro_tri);
G=zeros(2,2,nr_macro_tri);

% get all the transformations of the macrotriangles to the
% standard-triangle:
for k=1:nr_macro_tri
    tria_pts_x = pmacro_x(tmacro(k,:),:);
    tria_pts_y = pmacro_y(tmacro(k,:),:);
    [C(:,k), G(:,:,k)] = aff_trafo_glob2loc(tria_pts_x, tria_pts_y);
end

%check which point is in which macrotriangle
for macro_element_id = 1:nr_macro_tri
    pts_in_mac_tri = zeros(size(glob,1),1);
    C_big = repmat(C(:,macro_element_id),1,size(glob,1));
    pts_ref_dom = C_big + G(:,:,macro_element_id)*[pmic_x'; pmic_y']; %transform all points
    %check, which of the transformed points are in the reference triangle:
    % CAREFUL!!! See the "eps"! The comparison, wheter a value ist bigger 0
    % or smaller 1 is a bit sloppy...
    i=(pts_ref_dom(1,:)>=0-10*eps & pts_ref_dom(2,:)>=0-10*eps & pts_ref_dom(1,:)+pts_ref_dom(2,:)<=1+10*eps);
    %original (mathematically correct, but unfortunatelly not working
    %correctly...):
    %i=(pts_ref_dom(1,:)>=0 & pts_ref_dom(2,:)>=0 & pts_ref_dom(1,:)+pts_ref_dom(2,:)<=1);
    pts_in_mac_tri(i)=1; %index of the points in macro-triangle k
    
    % sort out those points, which are already in another triangle,
    % e.g. those which are on the edge of two triangles
    pts_in_mac_tri(pmic2tmac > 0) = 0;
    pmic2tmac = pmic2tmac + macro_element_id.* pts_in_mac_tri; 
end


if length(find(pmic2tmac == 0)) >= 1
    warning('for some points, no corresponding macrotriangle was found. maybe they are not in the domain')
end