function mp= Build_Cylinders()
% Bulid mp obj （contains models and all arguments)
    %% init
    pace=3.33;
    Pst=[-40 -40 0]-pace*0.52;
    Ped=[40 100 120];
    Ped=FCM_dV.round_Ped(Pst,Ped,pace);
    maxdeep=4;
    order=1;
    viewxy=[130 15];

    Shape.accessThreshold(pace/maxdeep/2/2/2);%
    Ped=FCM_dV.round_Ped(Pst,Ped,pace); 

    Oa=[0 0 0];
    Ra=[0 0 0];

    rvec=(Ped-Pst)/pace
    num_cells=rvec(1)*rvec(2)*rvec(3);% show num of cells
    fprintf("num_cells: %.0f\n",num_cells);% show num of cells

    %% add SS lib & python path
    additon_path=pwd;
    
    if count(py.sys.path,additon_path) == 0
        insert(py.sys.path,int32(0),additon_path);
    end

    %% assemble c0
    
    dh=20;
    dy=20;
    cy01=CSG('Cylinder',Oa,Ra,40,5,'cy01');
    cy02=CSG('Cylinder',Oa,Ra,27,39,'cy02');
    cy03=CSG('Cylinder',[0 0 34],Ra,35,56+dh,'cy03');
    cy04=CSG('Cylinder',[0 0 60],[-pi/2 0 0],25,70-dy,'cy04');

    ncy01=CSG('Cylinder',Oa,Ra,22.5,39,'ncy01');
    ncy02=CSG('Cylinder',[0 0 39],Ra,30,51+dh+4,'ncy02');
    ncy03=CSG('Cylinder',[0 0 60],[-pi/2 0 0],21,70,'ncy04');
    
    solidField=cy01+cy02+cy03+cy04;
    solidField.saveAs("solidField");
    emptyField=ncy01+ncy02+ncy03;
    emptyField.saveAs("emptyField");
    c=solidField-emptyField;


    V0=1;

    tic
    fcmv=FCM_V(Pst,Ped,pace,maxdeep,order,@c.calSDF,@(X) ones(size(X,1),1));% |||1 dV
    V0=fcmv.calIntegral_V();
    fprintf('V0: %f   ,Time: %f\n',V0,toc);

    %% build Dirichlet/Neumann boundary
    eepsD=2e-1;
    eepsN=1.5e-1;

    Dirich_csg=CSG('Cylinder',Oa,Ra,50,eps,'dir');
    Dirich_csg.judge_func=@(X) abs(X(:,3)-0)<eepsD;


    neu1=CSG('Cylinder',[0 0 90+dh],Ra,100,eps);
    neu1.saveAs("neu1");

    neu2=CSG('Cylinder',[0 50 60],[pi/2 0 0],50,eps);
    neu2.saveAs("neu2");
    neu1.judge_func=@(X) ((X(:,2)<40) & abs(X(:,3)-110)<eepsN);
    Neuman_csg=neu1+neu2;
    Neuman_csg.judge_func=@(X) ((X(:,2)<40) & abs(X(:,3)-110)<eepsN) | ((X(:,2)>=40) & abs(X(:,2)-50)<eepsN);
    Neuman_csg.fdv=@(X) (X(:,2)<40).*[0 0 -2e8*1e-3/1021.02]+(X(:,2)>=40).*[0 -1e8*1e-3/578.05 0];   
    Neuman_csg.eepsN=eepsN;

    fcmdv=FCM_dV(Pst,Ped,pace,maxdeep,order,@c.calSDF,@(X) ones(size(X,1),1));% |||1 dV
    S0=fcmdv.calIntegral_S(Dirich_csg)
    Sstd=3436.12;
    relaS=(S0-Sstd)/Sstd
    
    S1=fcmdv.calIntegral_S(Neuman_csg)
    S1std=1021.02+578.05;
    relaS1=(S1-S1std)/S1std
    
    %% show model
    clf;
    c=c+(CSG("NULL")*Dirich_csg);
    c.set_display_args(Pst,Ped,pace/2,'viewxy',viewxy);
    c.saveAs('c');
%     c.display_csg();
%     c.csg2stl('TMP_stp1');

    %% init mp(save handle)
    mp=ModelPlus(Pst,Ped,pace,order,maxdeep);
    mp.c=c;
    mp.Dirich_csg=Dirich_csg;
    mp.Neuman_csg=Neuman_csg;
    mp.viewxy=viewxy;
    mp.V0=V0;
    mp.E=2e11*1e-3;
    mp.u=0.3;

    
    if 1 && order<3
        saveKeH=SaveKe(mp,@compute_B_in_local,mp.pace,mp.order,mp);
        mp.saveKe=saveKeH;
    end

    %% set args here
    mp.fixedDir=1;
    mp.fixedNeu=0;
    mp.iscontrast=0;
    mp.isrecord=1;  

    mp.handleMap('neutop')=neu1;
    mp.handleMap('neuright')=neu2;
    %% init fmt
    tic
    fmt=FeatureMappingTable(c,Pst,Ped,pace);
    fmt.reflash_table();
    tini=toc;
    mp.recorder.time_maintain_fmt(end+1)=tini;
    mp.fmt=fmt;
    fprintf('num_node: %f\n',mp.id_cnt);


end