function [ sim,part ] = Initialize_mesh(sim,part,bound,corner,corner_num,center,partnum)
%INITIALIZE_MESH Summary of this function goes here
% XXX_d means discrete
%   Detailed explanation goes here
bound([1,2]) = bound([1,2])-sim.domain(1);
bound([3,4]) = bound([3,4])-sim.domain(3);
%in case of out of domain
if bound(1)<sim.domain(1)
    bound(1)=sim.domain(1);
end
if bound(3)<sim.domain(3)
    bound(3)=sim.domain(3);
end
if bound(2)>sim.domain(2)
    bound(2)=sim.domain(2);
end
if bound(4)>sim.domain(4)
    bound(4)=sim.domain(4);
end

bound_d = bound./sim.meshsize+1;
bound_d([1,3]) = floor(bound_d([1,3]));
bound_d([2,4]) = ceil(bound_d([2,4]));
bound_d(3) = sim.mesh_y_num+1-bound_d(3);
bound_d(4) = sim.mesh_y_num+1-bound_d(4);

sim.bound_d(partnum,:) = bound_d;

x = bound_d(2)-bound_d(1)+1;
y = bound_d(3)-bound_d(4)+1;
mesh = ones(y,x);
mesh_boundary = zeros(y,x);
mesh_calculate = zeros(y,x);
mesh_line_rec = zeros(y,x);
mesh_value = zeros(y,x,2);
mesh_value(:,:,:) = sim.mesh_value(bound_d(4):bound_d(3),...
    bound_d(1):bound_d(2),1:2);
boundary_rec = zeros(y,x,3);
boundary_rec(:,:,1:2) = mesh_value;

for xi=1:x
    change = 0;
    inter_line = 0;
    for yi=1:y
        for ii=1:corner_num
        inside = zeros(1,corner_num);
            point = [mesh_value(yi,xi,1);mesh_value(yi,xi,2)];
            if_inter = intersection...
                (corner(:,ii),corner(:,ii+1),center,point);
            inside(1,ii) = if_inter;
            if if_inter == 1
                inter_line = ii;
                break
            end
        end
        mesh_line_rec(yi,xi) = inter_line;
        if if_inter==0
            change = 1;
            break
        end
        mesh(yi,xi)=0;
    end
    
    if change==1 && yi~=y
        inter_line = 0;
        for yi=y:-1:1
            inside = zeros(1,corner_num);
            for ii=1:corner_num
                point = [mesh_value(yi,xi,1);mesh_value(yi,xi,2)];
                if_inter = intersection...
                    (corner(:,ii),corner(:,ii+1),center,point);
                inside(1,ii) = if_inter;
                if if_inter == 1
                    inter_line = ii;
                    break
                end
            end
            mesh_line_rec(yi,xi) = inter_line;
            if if_inter==0
                break
            end
            mesh(yi,xi)=0;
        end
    end
end

sim.mesh(bound_d(4):bound_d(3),bound_d(1):bound_d(2))=...
    sim.mesh(bound_d(4):bound_d(3),bound_d(1):bound_d(2))|mesh;

% get solid fraction
mesh_fraction = mesh; %0:air 0.2:boundary air 0.8:boundary solid 1:solid
for xi=1:x
    change = 0;
    for yi=1:y
        if mesh(yi,xi) == 1
            change = 1;
            point = [mesh_value(yi,xi,1);mesh_value(yi,xi,2)];
            mesh_boundary(yi,xi) = 1;
            %%%% solid fraction calculation part
            near_line = mesh_line_rec(yi,xi);
            mesh_fraction(yi,xi) = solid_fraction(sim,part,near_line,point,partnum);
            %%%%
            mesh_calculate(yi,xi) = 1;
            if yi>=2
                mesh_line_rec(yi-1,xi) = mesh_line_rec(yi,xi);
                %mesh_fraction(yi-1,xi) = 0;%calculate!!
                %%%%void solid fraction calculation part
                near_line = mesh_line_rec(yi-1,xi);
                point = [mesh_value(yi-1,xi,1);mesh_value(yi-1,xi,2)];
                mesh_fraction(yi-1,xi) = solid_fraction_void(sim,part,near_line,point,partnum);
                %%%%
                mesh_calculate(yi-1,xi) = 1;
            end
            break
        end
    end
    if change==1 && yi~=y
        for yi=y:-1:1
            if mesh(yi,xi)==1
                point = [mesh_value(yi,xi,1);mesh_value(yi,xi,2)];
                if mesh_boundary(yi,xi)==0
                    mesh_boundary(yi,xi) = 1;
                end
                mesh_line_rec(yi+1,xi) = mesh_line_rec(yi,xi);
                if mesh_calculate(yi,xi)==0
                    %%%% solid fraction calculation part
                    near_line = mesh_line_rec(yi,xi);
                    mesh_fraction(yi,xi) = solid_fraction(sim,part,near_line,point,partnum);
                    %%%%
                    mesh_calculate(yi,xi)=1;
                end
                if yi<y
                    %mesh_fraction(yi+1,xi)=0;%calculate!!
                    %%%%void solid fraction calculation part
                    near_line = mesh_line_rec(yi+1,xi);
                    point = [mesh_value(yi+1,xi,1);mesh_value(yi+1,xi,2)];
                    mesh_fraction(yi+1,xi) = solid_fraction_void(sim,part,near_line,point,partnum);
                    %%%%
                    mesh_calculate(yi+1,xi)=1;
                end
                break
            end
        end
    end
end
for yi=1:y
    change = 0;
    for xi=1:x
        if mesh(yi,xi)==1
            change = 1;
            point = [mesh_value(yi,xi,1);mesh_value(yi,xi,2)];
            if mesh_boundary(yi,xi)==0
                mesh_boundary(yi,xi) = 1;
            end
            if mesh_calculate(yi,xi)==0
                %%%% solid fraction calculation part
                near_line = 0;
                mesh_fraction(yi,xi) = solid_fraction(sim,part,near_line,point,partnum);
                %%%%
                mesh_calculate(yi,xi)=1;
            end
            
            if xi>=2
                mesh_line_rec(yi,xi-1) = mesh_line_rec(yi,xi);
                if mesh_calculate(yi,xi-1)==0
                    %mesh_fraction(yi,xi-1)=0;%calculate!!
                    %%%%void solid fraction calculation part
                    near_line = 0;
                    point = [mesh_value(yi,xi-1,1);mesh_value(yi,xi-1,2)];
                    mesh_fraction(yi,xi-1) = solid_fraction_void(sim,part,near_line,point,partnum);
                    %%%%
                    mesh_calculate(yi,xi-1)=1;
                end
            end
            break
        end
    end
    if change==1 && xi~=x
        for xi=x:-1:1
            if mesh(yi,xi)==1
                point = [mesh_value(yi,xi,1);mesh_value(yi,xi,2)];
                if mesh_boundary(yi,xi)==0
                    mesh_boundary(yi,xi) = 1;
                end
                mesh_line_rec(yi,xi+1) = mesh_line_rec(yi,xi);
                if mesh_calculate(yi,xi)==0
                    %%%% solid fraction calculation part
                    near_line = 0;
                    mesh_fraction(yi,xi) = solid_fraction(sim,part,near_line,point,partnum);
                    %%%%
                    mesh_calculate(yi,xi)=1;
                end
                if xi<x&&mesh_calculate(yi,xi+1)==0
                    %mesh_fraction(yi,xi+1)=0;%calculate!!
                    %%%%void solid fraction calculation part
                    near_line = 0;
                    point = [mesh_value(yi,xi+1,1);mesh_value(yi,xi+1,2)];
                    mesh_fraction(yi,xi+1) = solid_fraction_void(sim,part,near_line,point,partnum);
                    %%%%
                    mesh_calculate(yi,xi+1)=1;
                end
                break
            end
        end
    end
end
%save each particel's mesh information
mesh_fraction(mesh_fraction>(1-1e-10))=1;
mesh_fraction(mesh_fraction<1e-10)=0;
%%
fr = zeros(y,x);
fr(mesh_fraction>0&mesh_fraction<1)=1;
frr = zeros(y,x);
frr(mesh_fraction>0)=1;
boundary_rec(:,:,3) = fr;
name = ['part',num2str(partnum)];
part.mesh_boundary.(name) = boundary_rec;%%%

sim.mesh_solid_fra(bound_d(4):bound_d(3),bound_d(1):bound_d(2))=...
    sim.mesh_solid_fra(bound_d(4):bound_d(3),bound_d(1):bound_d(2))+mesh_fraction;

%%%%%%%%%%%%%%% update_velocity_mesh  sim.mesh_v
v_part = part.vel(:,partnum);
ang_vel = part.ang_vel(partnum);
mesh_vx = sim.mesh_v(bound_d(4):bound_d(3),bound_d(1):bound_d(2),1);
mesh_vy = sim.mesh_v(bound_d(4):bound_d(3),bound_d(1):bound_d(2),2);
v_calculated = sim.mesh_v_calculated(bound_d(4):bound_d(3),bound_d(1):bound_d(2),1);
for fri=1:y
    for frii=1:x
        if frr(fri,frii)==1
            poi = [mesh_value(fri,frii,1);mesh_value(fri,frii,2)]-center;
            poii = [-poi(2);poi(1)].*ang_vel;
            poi = poii+v_part;
            if v_calculated(fri,frii)==1
                mesh_vx(fri,frii) = (poi(1)+mesh_vx(fri,frii))/2;
                mesh_vy(fri,frii) = (poi(2)+mesh_vy(fri,frii))/2;
            else
                mesh_vx(fri,frii) = poi(1);
                mesh_vy(fri,frii) = poi(2);
            end
            v_calculated(fri,frii)=1;
        end
    end
end
sim.mesh_v(bound_d(4):bound_d(3),bound_d(1):bound_d(2),1) = mesh_vx;
sim.mesh_v(bound_d(4):bound_d(3),bound_d(1):bound_d(2),2) = mesh_vy;
sim.mesh_v_calculated(bound_d(4):bound_d(3),bound_d(1):bound_d(2),1) = v_calculated;
end