clc; clear; close all;

% zakladni promenne pro paralelizaci:
% partitions - cell array typu n_parts x 1 - prvky jsou vektory obsahuji
% globalni cisla elementu spadajicich do dane casti; pomoci teto promenne
% se rozhodne, ktere elementy patri na ktery procesor
% mapM_local, vmapM_local, mapP_local, vmapP_local, mapB_local,
% vmapB_local, mapM_interfaces_local, vmapM_interfaces_local,
% mapP_interfaces_local, vmapP_interfaces_local,
% neighbouring_partitions_local - vsechny cell array typu n_parts x 1 -
% jednolive prvky jsou pak prislusne maps, ktere se pouzivaji ve funkci
% ODEs_RHS a jou vytvoreny pro vsechny casti vypocetni site; dale posledni
% 3 promenne slouzi k volani dat z jineho PC z duvodu, ze prilehajici body
% k bodum na hranicich casti site lezi jiz v jine casti a jsou tedy
% uchovavany na jinem PC (pricemz prislusne mapP a vmapP jiz primo obsahuje
% pristupove informace k potrebnym datum na jinem PC) 

% asi preznacit mapM na mapM_global atd.

visual_debug = 1;
debug_part_index = 1;

Globals3D_modified;

N = 1;

[P,E,E_parts,E_parts_neighs] = MshToMatlab3D_with_partitions('Advection3D_h=2.00_2partitions.msh');

StartUp3D_modified;

n_parts = max(E_parts);

partitions_matrix = zeros(n_parts,K);
pointers = ones(1,n_parts);
for i = 1:K
    part_ID = E_parts(i);
    partitions_matrix(part_ID,pointers(part_ID)) = i;
    pointers(part_ID) = pointers(part_ID) + 1;
end
partitions = cell(n_parts,1);
for i = 1:n_parts
    partitions_vector = partitions_matrix(i,:);
    partitions_vector(pointers(i):end) = [];
    partitions{i} = partitions_vector;
end

% mapM, vmapM, mapP, vmapP - v globalni siti
%mapM_local = cell(n_parts,1);
%vmapM_local = cell(n_parts,1);
%mapP_local = cell(n_parts,1);
%vmapP_local = cell(n_parts,1);
%for i = 1:n_parts
%    elems_in_partition = partitions{i};
%    K_local = length(elems_in_partition);
%    n_local = K_local * N_faces * N_fp;
%    mapM_local_temp = zeros(n_local,1);
%    vmapM_local_temp = zeros(n_local,1);
%    mapP_local_temp = zeros(n_local,1);
%    vmapP_local_temp = zeros(n_local,1);
%    for j = 1:K_local
%        indeces = ((j - 1) * N_faces * N_fp + 1):(j * N_faces * N_fp);
%        indeces_global = ((elems_in_partition(j) - 1) * N_faces * N_fp + 1):(elems_in_partition(j) * N_faces * N_fp);
%        mapM_local_temp(indeces) = mapM(indeces_global);
%        vmapM_local_temp(indeces) = vmapM(indeces_global);
%        mapP_local_temp(indeces) = mapP(indeces_global);
%        vmapP_local_temp(indeces) = vmapP(indeces_global);
%    end
%    mapM_local{i} = mapM_local_temp;
%    vmapM_local{i} = vmapM_local_temp;
%    mapP_local{i} = mapP_local_temp;
%    vmapP_local{i} = vmapP_local_temp;
%end
%if long_maps_debug
%temp = mapM_local{debug_part_index};
%temp'
%temp = mapP_local{debug_part_index};
%temp'
%end
%% mapM, vmapM, mapP, vmapP - v lokalni siti
%for i = 1:n_parts
%    elems_in_partition = partitions{i};
%    mapM_partition = mapM_local{i};
%    vmapM_partition = vmapM_local{i};
%    mapP_partition = mapP_local{i};
%    vmapP_partition = vmapP_local{i};
%    n_local = length(mapM_partition);
%    mapM_local_temp = zeros(n_local,1);
%    vmapM_local_temp = zeros(n_local,1);
%    mapP_local_temp = zeros(n_local,1);
%    vmapP_local_temp = zeros(n_local,1);
%    for j = 1:n_local
%        mapM_global_ID = mapM_partition(j);
%        mapM_local_ID = GlobalToLocal_map(mapM_global_ID,elems_in_partition);
%        vmapM_global_ID = vmapM_partition(j);
%        vmapM_local_ID = GlobalToLocal_vmap(vmapM_global_ID,elems_in_partition);
%        mapP_global_ID = mapP_partition(j);
%        mapP_local_ID = GlobalToLocal_map(mapP_global_ID,elems_in_partition);
%        vmapP_global_ID = vmapP_partition(j);
%        vmapP_local_ID = GlobalToLocal_vmap(vmapP_global_ID,elems_in_partition);
%        
%        mapM_local_temp(j) = mapM_local_ID;
%        vmapM_local_temp(j) = vmapM_local_ID;
%        if isempty(mapP_local_ID)
%            mapP_local_temp(j) = -mapM_local_ID;
%        else
%            mapP_local_temp(j) = mapP_local_ID;
%        end
%        if isempty(vmapP_local_ID)
%            vmapP_local_temp(j) = -vmapM_local_ID;
%        else
%            vmapP_local_temp(j) = vmapP_local_ID;
%        end
%    end
%    mapM_local{i} = mapM_local_temp;
%    vmapM_local{i} = vmapM_local_temp;
%    mapP_local{i} = mapP_local_temp;
%    vmapP_local{i} = vmapP_local_temp;
%end
%if long_maps_debug
%temp = mapM_local{debug_part_index};
%temp'
%temp = mapP_local{debug_part_index};
%temp'
%end
%
%% urceni bodu, ktere lezi na fyzikalni hranici
%mapB_local = cell(n_parts,1);
%vmapB_local = cell(n_parts,1);
%for i = 1:n_parts
%    mapB_local{i} = mapM_local{i}(mapM_local{i} == mapP_local{i});
%    vmapB_local{i} = vmapM_local{i}(mapM_local{i} == mapP_local{i});
%end
%
%% urceni bodu, ktere lezi uvnitr dane partition (nesousedi s jinou partition)
%mapM_inner_local = cell(n_parts,1);
%vmapM_inner_local = cell(n_parts,1);
%mapP_inner_local = cell(n_parts,1);
%vmapP_inner_local = cell(n_parts,1);
%for i = 1:n_parts
%    mapM_inner_local{i} = mapM_local{i}(mapM_local{i} ~= -mapP_local{i});
%    vmapM_inner_local{i} = vmapM_local{i}(mapM_local{i} ~= -mapP_local{i});
%    mapP_inner_local{i} = mapP_local{i}(mapM_local{i} ~= -mapP_local{i});
%    vmapP_inner_local{i} = vmapP_local{i}(mapM_local{i} ~= -mapP_local{i});
%end
%if long_maps_debug
%temp = mapM_inner_local{debug_part_index};
%temp'
%temp = mapP_inner_local{debug_part_index};
%temp'
%end
%
%% urceni bodu na hranicich mezi castmi site
%mapM_interfaces_local = cell(n_parts,1);
%vmapM_interfaces_local = cell(n_parts,1);
%for i = 1:n_parts
%    mapM_interfaces_local{i} = mapM_local{i}(mapM_local{i} == -mapP_local{i});
%    vmapM_interfaces_local{i} = vmapM_local{i}(mapM_local{i} == -mapP_local{i});
%end
%if maps_debug
%temp = vmapM_interfaces_local{debug_part_index};
%temp'
%end
%neighbouring_partitions_IDs_local = cell(n_parts,1);
%for i = 1:n_parts
%    elems_in_partition = partitions{i};
%    mapM_interfaces_partition = mapM_interfaces_local{i};
%    n_local = length(mapM_interfaces_partition);
%    neigh_parts_temp = zeros(n_local,1);
%    for j = 1:n_local        
%        mapM_local_ID = mapM_interfaces_partition(j);
%        [mapM_global_ID,elem_global,neigh_face] = LocalToGlobal_map(mapM_local_ID,elems_in_partition);
%        
%        neigh_elem_global = EtoE(elem_global,neigh_face);
%        neigh_part_temp = E_parts(neigh_elem_global);
%       
%        neigh_parts_temp(j) = neigh_part_temp;
%    end
%    neighbouring_partitions_IDs_local{i} = neigh_parts_temp;
%end
%if maps_debug
%temp = neighbouring_partitions_IDs_local{debug_part_index};
%temp'
%end
%neighbouring_partitions_IDs_unique_local = cell(n_parts,1);
%for i = 1:n_parts
%    neighbouring_partitions_IDs_unique_local{i} = unique(neighbouring_partitions_IDs_local{i});
%end
%temp = neighbouring_partitions_IDs_unique_local{debug_part_index};
%temp'
%% spocitani permutacnich vektoru, aby mapM_interfaces a k nim
%% predani sousede byly pak na sebe nalehajici body
%permutation_vector_local = cell(n_parts,1);
%for i = 1:n_parts
%    elems_in_partition = partitions{i};
%    mapM_interfaces_partition = mapM_interfaces_local{i};
%    n_local = length(mapM_interfaces_partition);
%    mapM_interfaces_global = zeros(n_local,1);
%    for j = 1:n_local
%        mapM_local_ID = mapM_interfaces_partition(j);
%        [mapM_global_ID,elem_global,neigh_face] = LocalToGlobal_map(mapM_local_ID,elems_in_partition);
%        mapM_interfaces_global(j) = mapM_global_ID;
%    end
%    mapM_interface_global_IDs = mapP(mapM_interfaces_global);
%    global_neighbours = zeros(n_local,1);
%    current_neighbouring_partitions = neighbouring_partitions_IDs_local{i};
%    current_neighbouring_partitions_unique = neighbouring_partitions_IDs_unique_local{i};
%    n_neighbouring_partitions = length(current_neighbouring_partitions_unique);
%    for j = 1:n_neighbouring_partitions
%        partition_ID = current_neighbouring_partitions_unique(j);
%        elems_in_neighbouring_partition = partitions{partition_ID};
%        indeces = find(current_neighbouring_partitions == partition_ID);
%        mapM_interfaces_neighbouring_partition = mapM_interfaces_local{partition_ID}(neighbouring_partitions_IDs_local{partition_ID} == i);
%        n2_local = length(mapM_interfaces_neighbouring_partition);
%        mapM_interfaces_neighbouring_partition_global = zeros(n2_local,1);
%        for k = 1:n2_local
%            mapM_local_ID = mapM_interfaces_neighbouring_partition(k);
%            [mapM_global_ID,elem_global,neigh_face] = LocalToGlobal_map(mapM_local_ID,elems_in_neighbouring_partition);
%            mapM_interfaces_neighbouring_partition_global(k) = mapM_global_ID;
%        end
%        global_neighbours(indeces) = mapM_interfaces_neighbouring_partition_global;
%    end
%    permutation_vector_temp = zeros(n_local,1);
%    for j = 1:n_local
%        temp_ID = mapM_interface_global_IDs(j);
%        position = find(global_neighbours == temp_ID);
%        permutation_vector_temp(j) = position;
%    end
%    permutation_vector_local{i} = permutation_vector_temp;
%end



% vykresleni vsech elementu barevne dle casti site, do ktere patri
if visual_debug == 1
    figure;
    heights = linspace(-1,1,n_parts);
    for i = 1:n_parts
        elements = partitions{i};
        E_temp = E(:,elements);
        tetramesh(E_temp',P',heights(i) * ones(1,length(elements)));
    end
    axis equal;
    xlim([min(P(1,:)), max(P(1,:))]); ylim([min(P(2,:)), max(P(2,:))]); zlim([min(P(3,:)), max(P(3,:))]);
    title('casti site z Gmsh');
end

%% test mapM_local, vmapM_local
%if visual_debug == 1
%    elems_in_partition = partitions{debug_part_index};
%    x_partition = x(:,partitions{debug_part_index});
%    y_partition = y(:,partitions{debug_part_index});
%    Fx_partition = Fx(:,partitions{debug_part_index});
%    Fy_partition = Fy(:,partitions{debug_part_index});
%    marker_type = '*';
%    marker_size = 12;
%    figure; % test vmapM_local
%    plot(x_partition(vmapM_local{debug_part_index}),y_partition(vmapM_local{debug_part_index}),marker_type,'MarkerSize',marker_size); hold on;
%    trimesh(E(:,elems_in_partition)',P(1,:),P(2,:),'Color','black');
%    axis equal;
%    title('test vmapM local');
%    xlim([min(P(1,:)), max(P(1,:))]); ylim([min(P(2,:)), max(P(2,:))]);
%    figure; % test mapM_local
%    plot(Fx_partition(mapM_local{debug_part_index}),Fy_partition(mapM_local{debug_part_index}),marker_type,'MarkerSize',marker_size); hold on;
%    trimesh(E(:,elems_in_partition)',P(1,:),P(2,:),'Color','black');
%    axis equal;
%    xlim([min(P(1,:)), max(P(1,:))]); ylim([min(P(2,:)), max(P(2,:))]);
%    title('test mapM local');
%end
%
%% test mapB_local, vmapB_local
%if visual_debug == 1
%    elems_in_partition = partitions{debug_part_index};
%    x_partition = x(:,partitions{debug_part_index});
%    y_partition = y(:,partitions{debug_part_index});
%    Fx_partition = Fx(:,partitions{debug_part_index});
%    Fy_partition = Fy(:,partitions{debug_part_index});
%    marker_type = '*';
%    marker_size = 12;
%    figure; % test vmapB_local
%    plot(x_partition(vmapB_local{debug_part_index}),y_partition(vmapB_local{debug_part_index}),marker_type,'MarkerSize',marker_size); hold on;
%    trimesh(E(:,elems_in_partition)',P(1,:),P(2,:),'Color','black');
%    axis equal;
%    xlim([min(P(1,:)), max(P(1,:))]); ylim([min(P(2,:)), max(P(2,:))]);
%    title('test vmapB local');
%    figure; % test mapB_local
%    plot(Fx_partition(mapB_local{debug_part_index}),Fy_partition(mapB_local{debug_part_index}),marker_type,'MarkerSize',marker_size); hold on;
%    trimesh(E(:,elems_in_partition)',P(1,:),P(2,:),'Color','black');
%    axis equal;
%    xlim([min(P(1,:)), max(P(1,:))]); ylim([min(P(2,:)), max(P(2,:))]);
%    title('test mapB local');
%end
%
%% test mapM_inner_local, vmapM_inner_local
%if visual_debug == 1
%    elems_in_partition = partitions{debug_part_index};
%    x_partition = x(:,partitions{debug_part_index});
%    y_partition = y(:,partitions{debug_part_index});
%    Fx_partition = Fx(:,partitions{debug_part_index});
%    Fy_partition = Fy(:,partitions{debug_part_index});
%    marker_type = '*';
%    marker_size = 12;
%    figure; % test vmapB_local
%    plot(x_partition(vmapM_inner_local{debug_part_index}),y_partition(vmapM_inner_local{debug_part_index}),marker_type,'MarkerSize',marker_size); hold on;
%    trimesh(E(:,elems_in_partition)',P(1,:),P(2,:),'Color','black');
%    axis equal;
%    xlim([min(P(1,:)), max(P(1,:))]); ylim([min(P(2,:)), max(P(2,:))]);
%    title('test vmapM inner local');
%    figure; % test mapB_local
%    plot(Fx_partition(mapM_inner_local{debug_part_index}),Fy_partition(mapM_inner_local{debug_part_index}),marker_type,'MarkerSize',marker_size); hold on;
%    trimesh(E(:,elems_in_partition)',P(1,:),P(2,:),'Color','black');
%    axis equal;
%    xlim([min(P(1,:)), max(P(1,:))]); ylim([min(P(2,:)), max(P(2,:))]);
%    title('test mapM inner local');
%end
%
%% test mapM_interfaces_local, vmapM_interfaces_local a neighbouring_partitions_IDs_local
%if visual_debug == 1
%    
%    elems_in_partition = partitions{debug_part_index};
%    x_partition = x(:,partitions{debug_part_index});
%    y_partition = y(:,partitions{debug_part_index});
%    Fx_partition = Fx(:,partitions{debug_part_index});
%    Fy_partition = Fy(:,partitions{debug_part_index});
%    marker_type = '*';
%    marker_type2 = 'o';
%    marker_size = 12;
%    
%    neigh_parts = neighbouring_partitions_IDs_local{debug_part_index};
%    n_neigh_parts = length(neighbouring_partitions_IDs_unique_local{debug_part_index});
%    
%    figure; % test vmapM_interfaces_local
%    plot(x_partition(vmapM_interfaces_local{debug_part_index}),y_partition(vmapM_interfaces_local{debug_part_index}),marker_type,'MarkerSize',marker_size); hold on;
%    for i = 1:n_neigh_parts
%        neigh_part_ID = neighbouring_partitions_IDs_unique_local{debug_part_index}(i);
%        neigh_parts_on_neigh_part = neighbouring_partitions_IDs_local{neigh_part_ID};
%        x_neigh_part = x(:,partitions{neigh_part_ID});
%        y_neigh_part = y(:,partitions{neigh_part_ID});
%        plot(x_neigh_part(vmapM_interfaces_local{neigh_part_ID}(neigh_parts_on_neigh_part == debug_part_index)),y_neigh_part(vmapM_interfaces_local{neigh_part_ID}(neigh_parts_on_neigh_part == debug_part_index)),marker_type2,'MarkerSize',marker_size); hold on;
%    end
%    trimesh(E(:,elems_in_partition)',P(1,:),P(2,:),'Color','black');
%    axis equal;
%    xlim([min(P(1,:)), max(P(1,:))]); ylim([min(P(2,:)), max(P(2,:))]);
%    title('test vmapM and vmapP interfaces local');
%    figure; % test mapM_interfaces_local
%    plot(Fx_partition(mapM_interfaces_local{debug_part_index}),Fy_partition(mapM_interfaces_local{debug_part_index}),marker_type,'MarkerSize',marker_size); hold on;
%    for i = 1:n_neigh_parts
%        neigh_part_ID = neighbouring_partitions_IDs_unique_local{debug_part_index}(i);
%        neigh_parts_on_neigh_part = neighbouring_partitions_IDs_local{neigh_part_ID};
%        Fx_neigh_part = Fx(:,partitions{neigh_part_ID});
%        Fy_neigh_part = Fy(:,partitions{neigh_part_ID});
%        plot(Fx_neigh_part(mapM_interfaces_local{neigh_part_ID}(neigh_parts_on_neigh_part == debug_part_index)),Fy_neigh_part(mapM_interfaces_local{neigh_part_ID}(neigh_parts_on_neigh_part == debug_part_index)),marker_type2,'MarkerSize',marker_size); hold on;
%    end
%    trimesh(E(:,elems_in_partition)',P(1,:),P(2,:),'Color','black');
%    axis equal;
%    xlim([min(P(1,:)), max(P(1,:))]); ylim([min(P(2,:)), max(P(2,:))]);
%    title('test mapM and mapP interfaces local');
%end
%
%
%% debug jednolivych lokalnich vektoru a jejich delek
%fprintf(1,'mapM_local{%d} ma %d prvku\n',debug_part_index,length(mapM_local{debug_part_index}));
%fprintf(1,'mapM_inner_local{%d} ma %d prvku\n',debug_part_index,length(mapM_inner_local{debug_part_index}));
%fprintf(1,'mapM_interfaces_local{%d} ma %d prvku\n',debug_part_index,length(mapM_interfaces_local{debug_part_index}));
%fprintf(1,'mapM_interfaces_local{%d} = \n',debug_part_index);
%temp = mapM_interfaces_local{debug_part_index};
%temp'
%fprintf(1,'mapB_local{%d} ma %d prvku\n',debug_part_index,length(mapB_local{debug_part_index}));
%fprintf(1,'mapB_local{%d} = \n',debug_part_index);
%temp = mapB_local{debug_part_index};
%temp'




% NEPOTREBNY ZDROJAK

%IDs = ~cellfun('isempty',E_parts_neighs);
%E_parts_neighs_full = E_parts_neighs(IDs);
%% vykresleni tech elementu, ktere dle Gmsh sousedi s elementem v jine casti site
%if visual_debug == 1
%    figure;
%    heights = linspace(-1,1,n_parts);
%    for i = 1:n_parts
%        elements = partitions{i};
%        E_temp = E(:,elements);
%        IDs_temp = IDs(elements);
%        trisurf(E_temp(:,IDs_temp)',P(1,:),P(2,:),heights(i) * ones(1,N_v),'LineWidth',1.5); hold on;
%    end
%    view(0, 90);
%    axis equal;
%    xlim([min(P(1,:)), max(P(1,:))]); ylim([min(P(2,:)), max(P(2,:))]);
%    title('elementy, ktere lezi na hranicich mezi castmi site dle Gmsh');
%end
%
%
%% je videt, ze Gmsh chape jako souseda i element, ktery ma spolecny pouze vrchol
%% v DGM uvazujeme, ze sousede jsou elementy, ktere maji spolecnou hranu ve 2D
%E_neigh_parts_faces = cell(1,K);
%E_parts_neighs_modif = E_parts_neighs;
%neighs_found = 0;
%for i = find(IDs)
%    elem_def = EtoE(i,:);
%    neighs_DG = setdiff(elem_def,i);
%    parts_neighs = E_parts_neighs{i};
%    n_parts_neighs = length(parts_neighs);
%    parts_neighs_modif = parts_neighs;
%    neigh_parts_faces = parts_neighs;
%    for j = 1:n_parts_neighs
%        neight_part = parts_neighs(j);
%        elems_neigh_part = find(E_parts == neight_part);
%        common_elem = intersect(elems_neigh_part,neighs_DG);
%        if ~isempty(common_elem)
%            neighs_found = neighs_found + 1;
%            parts_neighs_modif(neighs_found) = neight_part;
%            neigh_parts_faces(neighs_found) = find(elem_def == common_elem);
%        end
%    end
%    if (neighs_found < n_parts_neighs)
%        parts_neighs_modif((neighs_found + 1):end) = [];
%        neigh_parts_faces((neighs_found + 1):end) = [];
%    end
%    E_parts_neighs_modif{i} = parts_neighs_modif;
%    E_neigh_parts_faces{i} = neigh_parts_faces;
%    
%    neighs_found = 0;
%end
%
%
%% vykresleni spravnych elementu, ktere jsou na hranicich casti site
%IDs_modif = ~cellfun('isempty',E_parts_neighs_modif);
%elem_IDs = find(IDs_modif);
%E_neigh_parts_faces_full_modif = E_neigh_parts_faces(IDs_modif);
%E_parts_neighs_full_modif = E_parts_neighs_modif(IDs_modif);
%if visual_debug == 1
%    figure;
%    heights = linspace(-1,1,n_parts);
%    for i = 1:n_parts
%        elements = partitions{i};
%        E_temp = E(:,elements);
%        IDs_temp = IDs_modif(elements);
%        trisurf(E_temp(:,IDs_temp)',P(1,:),P(2,:),heights(i) * ones(1,N_v),'LineWidth',1.5); hold on;
%    end
%    face_def = [1, 2; 2, 3; 1, 3];
%    for i = 1:length(elem_IDs)
%        elem_ID = elem_IDs(i);
%        faces = E_neigh_parts_faces_full_modif{i};
%        for j = 1:length(faces)
%            face = face_def(faces(j),:);
%            x_face = [P(1,E(face(1),elem_ID)),P(1,E(face(2),elem_ID))];
%            y_face = [P(2,E(face(1),elem_ID)),P(2,E(face(2),elem_ID))];
%            plot3(x_face,y_face,[2, 2],'LineWidth',3,'Color','red'); hold on;
%        end
%    end
%    view(0, 90);
%    axis equal;
%    xlim([min(P(1,:)), max(P(1,:))]); ylim([min(P(2,:)), max(P(2,:))]);
%    title('spravne elementy, ktere lezi na hranicich mezi castmi site');
%end
%
%% lokalni IDs elementu na hranici v kazde casti site
%boundary_elems_IDs_global = cell(n_parts,1);
%boundary_elems_IDs_local = cell(n_parts,1);
%for i = 1:n_parts
%    boundary_elems_IDs_global{i} = intersect(partitions{i},elem_IDs);
%    boundary_elems_IDs_local{i} = find(ismember(partitions{i},elem_IDs));
%end


% urceni mapB a vmapB - udelano efektivneji
%% mapB and vmapB - v globalni siti
%mapB_cell = cell(1,K);
%vmapB_cell = cell(1,K);
%for i = 1:K
%    mapB_cell{i} = intersect(((i - 1) * N_faces * N_fp + 1):(i * N_faces * N_fp),mapB);
%    vmapB_cell{i} = intersect(((i - 1) * N_p + 1):(i * N_p),vmapB);
%end
%mapB_local = cell(n_parts,1);
%vmapB_local = cell(n_parts,1);
%for i = 1:n_parts
%    K_local = length(partitions{i});
%    n_boundary_points = 0;
%    for j = 1:K_local
%        n_boundary_points = n_boundary_points + length(mapB_cell{partitions{i}(j)});
%    end
%    mapB_local_temp = zeros(n_boundary_points,1);
%    vmapB_local_temp = zeros(n_boundary_points,1);
%    pointer_end = 0;
%    for j = 1:K_local
%        n_boundary_points_local = length(mapB_cell{partitions{i}(j)});
%        if (n_boundary_points_local > 0)
%            mapB_local_temp((pointer_end + 1):(pointer_end + n_boundary_points_local)) = mapB_cell{partitions{i}(j)};
%            vmapB_local_temp((pointer_end + 1):(pointer_end + n_boundary_points_local)) = vmapB_cell{partitions{i}(j)};
%            pointer_end = pointer_end + n_boundary_points_local;
%        end
%    end
%    mapB_local{i} = mapB_local_temp;
%    vmapB_local{i} = vmapB_local_temp;
%end
%% mapB and vmapB - v lokalni siti
%for i = 1:n_parts
%    elems_in_partition = partitions{i};
%    mapB_partition = mapB_local{i};
%    vmapB_partition = vmapB_local{i};
%    n_local = length(mapB_local{i});
%    mapB_local_temp = zeros(n_local,1);
%    vmapB_local_temp = zeros(n_local,1);
%    for j = 1:n_local
%        mapB_global_ID = mapB_partition(j);
%        mapB_local_ID = GlobalToLocal_map(mapB_global_ID,elems_in_partition);
%        vmapB_global_ID = vmapB_partition(j);
%        vmapB_local_ID = GlobalToLocal_vmap(vmapB_global_ID,elems_in_partition);
%        
%        mapB_local_temp(j) = mapB_local_ID;
%        vmapB_local_temp(j) = vmapB_local_ID;
%    end
%    mapB_local{i} = mapB_local_temp;
%    vmapB_local{i} = vmapB_local_temp;
%end