function mp= Build_BearingBracket()
% Bulid mp obj （contains models and all arguments)
    %% init
    pace=2.5;
    Pst=[-6 -20 -15]-pace*0.52;
    Ped=[116 20 80];
    Ped=FCM_dV.round_Ped(Pst,Ped,pace); 
    maxdeep=4;
    order=2;
    viewxy=[-30 30];
    
    Shape.accessThreshold(pace/maxdeep/2/2/2);%
    Oa=[0 0 0];
    Ra=[0 0 0];

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


    %% assemble c0
    c=CSG("STL",[pwd,'/asset/SolidWoks Files/BearingBracket/BearingB.STL'],'init');

    num_c_in_cbs=24;
    cbs_cvs=zeros(num_c_in_cbs,1)+4;
%     cur_cbs=CBS(cbs_cvs,3,[0 0 0]);

    center_Y=[102.54 0 25.83
            101.27 0 42.08
            83.09 0 59.03
            71.3 0 48.73
            60.51 0 38.68
            32.14 0 9.35
%             30.39 0 7.35
        ];
    center_Z=[
        84.89 0 11.86 
        71.3 0 11.86
        56.51 0 11.86
        18.38 0 -5
        ];
    center_X=[
        41.83 0 17.55
        ];
    
    num_X=size(center_X,1);
    num_Y=size(center_Y,1);
    num_Z=size(center_Z,1);
    num_cbs=num_X+num_Y+num_Z;
%     skt_list=[];

    cbs_list=[];
    %---------------------------for Y
    heg_Y=40;
    for i_Y=1:num_Y
        cur_cbs=CBS(cbs_cvs,3,[0 0 0]);
        cbs_list=[cbs_list cur_cbs];
        skt_list(i_Y,1)=Sketch(center_Y(i_Y,:)-[0 heg_Y/2 0],[-pi/2 0 0],heg_Y,cur_cbs);
    end
    %---------------------------for Z
    heg_Z=100;
    for i_Z=1:size(center_Z,1)
        cur_cbs=CBS(cbs_cvs,3,[0 0 0]);
        cbs_list=[cbs_list cur_cbs];
        skt_list(num_Y+i_Z,1)=Sketch(center_Z(i_Z,:),[0 0 0],heg_Z,cur_cbs);
    end
    %---------------------------for X
    heg_X=45;
    for i_X=1:size(center_X,1)
        cur_cbs=CBS(cbs_cvs,3,[0 0 0]);
        cbs_list=[cbs_list cur_cbs];
        skt_list(num_Y+num_Z+i_X,1)=Sketch(center_X(i_X,:)-[1 0 0],[0 pi/2 0],heg_X,cur_cbs);
    end
    

    %---------------------------skt2csg
    skt_csg_list=[];
    for i_skt=1:length(skt_list)
        cur_csg=CSG();cur_csg.init(skt_list(i_skt));
        skt_csg_list=[skt_csg_list cur_csg];
    end
    csg_list=skt_csg_list;
    %---------------------------limit csg
    pla_csg_top=CSG();pla_csg_top.init(Plane3D([83.14 -12.91 47.37],[1 0 -1]));
    pla_csg_down=CSG();pla_csg_down.init(Plane3D([94.23 3.56 9.66],[0 0 1]));
    for i_csg=1:2
        cur_csg=csg_list(i_csg);
        cur_csg=cur_csg*pla_csg_top*pla_csg_down;
        csg_list(i_csg)=cur_csg;
    end


    pla_csg_left=CSG();pla_csg_left.init(Plane3D([93.23 10.95 42.77],[-1 0 0]));
    pla_csg_right=CSG();pla_csg_right.init(Plane3D([47.62 18.97 11.86],[1 0 0]));
    for i_csg=[3:5,7:9]
        cur_csg=csg_list(i_csg);
        cur_csg=cur_csg*pla_csg_left*pla_csg_right;
        csg_list(i_csg)=cur_csg;
    end       

    front_csg=csg_list(6);
    pla_csg_down=CSG();pla_csg_down.init(Plane3D([0 0 2.91],[0 0 1]));
    pla_csg_left=CSG();pla_csg_left.init(Plane3D([41.83 13.25 2.91],[-1 0 0]));
    pla_csg_3=CSG();pla_csg_3.init(Plane3D([45.23 19.05 30],[1 0 -1]));
    front_csg=front_csg*pla_csg_down*pla_csg_left*pla_csg_3;
    csg_list(6)=front_csg;
    
    right_csg=csg_list(10);
    cir1=Circle([36.83 -8.25 0],5);
    cir2=Circle([36.83 8.25 0],5);
    cir3=Circle([0 -9.52 0],5);
    cir4=Circle([0 9.52 0],5);
    limit_skt=Sketch([0 0 -5],[0 0 0],8,(cir1+cir2+cir3+cir4));
    limit_csg=CSG();limit_csg.init(limit_skt);

    pla_csg_down=CSG();pla_csg_down.init(Plane3D([0.08 14.52 2.91],[-0.5 -14.52 0]));
    pla_csg_up=CSG();pla_csg_up.init(Plane3D([0.08 -14.52 2.91],[-0.5 14.52 0]));
    pla_csg_left=CSG();pla_csg_left.init(Plane3D([36.83 8.26 2.91],[-1 0 0]));
    pla_csg_right=CSG();pla_csg_right.init(Plane3D([0 -9.52 2.91],[1 0 0]));    
    csg_list(10)=(right_csg-limit_csg)*pla_csg_down*pla_csg_up*pla_csg_left*pla_csg_right;

    last_csg=csg_list(end);
    pla_csg_down=CSG();pla_csg_down.init(Plane3D([41.83 8.26 2.91],[0 0 1]));
    pla_csg_up=CSG();pla_csg_up.init(Plane3D([41.83 8.26 34.35],[0 0 -1]));
    
    pla_csg_left=CSG();pla_csg_left.init(Plane3D([0.08 14.52 2.91],[-0.5 -14.52 0]));
    pla_csg_right=CSG();pla_csg_right.init(Plane3D([0.08 -14.52 2.91],[-0.5 14.52 0]));    
    last_csg=last_csg*pla_csg_down*pla_csg_up*pla_csg_left*pla_csg_right;
    csg_list(end)=last_csg;    
    
    %---------------------------assemble c
    sub_csgs=csg_list(1);
    for i_csg=2:length(csg_list)
        sub_csgs=sub_csgs+csg_list(i_csg);
    end
    
    
   
    c=c-sub_csgs;
    % show model
    clf;
    c.set_display_args([Pst(1) Pst(2) Pst(3)],Ped,pace/5,'viewxy',viewxy);
    c.display_csg();

    V0=1;
    % time and calculate Volume
    tic
    fcmv=FCM_V(Pst,Ped,pace,maxdeep,order,@c.calSF,@(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=1e-1;
eepsN=8e-1;

DC2d=Circle([36.83 -8.26 0],3)+Circle([36.83 8.26 0],3)+Circle([0 -9.52 0],3)+Circle([0 9.52 0],3);
Dskt=Sketch([0 0 -5],[0 0 0],5+6+5,DC2d);
Dirich_csg=CSG();Dirich_csg.init(Dskt);
Dirich_csg.saveAs('Dir');
Dirich_csg.judge_func=@(X)  X(:,1)<18 & abs(vecnorm([X(:,1),abs(X(:,2))]-[0 9.53],2,2)-3)<eepsD | X(:,1)>=18 & abs(vecnorm([X(:,1),abs(X(:,2))]-[36.83 8.26],2,2)-3)<eepsD;
%------------------------------------------Neu

Ncir=Circle([0 0 0],7.065);
Nskt=Sketch([101.6 0 -2.21]-[0 5 0],[-pi/2 0 0],10,Ncir);
Neuman_csg=CSG();
Neuman_csg.init(Nskt,'Neu');
Neuman_csg.judge_func=@(X) abs(vecnorm(X(:,[1,3])-[101.6 -2.21],2,2)-7.065)<eepsN;
Neuman_csg.fdv=@cal_NeuF;

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=114.64*4;
relaSD=(S0-Sstd)/Sstd

S1=fcmdv.calIntegral_S(Neuman_csg)
S1std=373.14;
relaSN=(S1-S1std)/S1std 

    %% show model
    clf;
    c.set_display_args([Pst(1) 0 Pst(3)],Ped,pace/2,'viewxy',viewxy);
    c.saveAs('c');
    c.display_csg();



    %% 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.V0=V0;
    mp.viewxy=viewxy;
    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=1;
    mp.iscontrast=0;
    mp.isrecord=1;
    mp.frac=0.6;
    mp.maxloop=100-100+50;
    mp.interp_eeps=pace/(2^(maxdeep-1));
%     mp.Neu_norm_flg=1;
    %% save handle here


    
    %% 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);

    %% init x
    % x y num_c_in_cbs-2: tot num_c_in_cbs * 10=100
    x0=zeros(num_c_in_cbs,1);
    lb=x0;
    ub=x0;
    for i_x=1:num_cbs
        st_idx=(i_x-1)*num_c_in_cbs;
        for ii_x=1:num_c_in_cbs
            if ii_x<=2
                % x y posision
                
                % init
                lb(st_idx+ii_x)=-20;
                ub(st_idx+ii_x)=20;

                % fix x y: 3,4,5, 7,8,9
                if i_x>=3 && i_x<=5 || i_x>=7 && i_x<=9
                    lb(st_idx+ii_x)=-1;
                    ub(st_idx+ii_x)=1;
                end

                
                % fix x: 1,2
                if i_x>=0 && i_x<=2
                    if ii_x==1
                        lb(st_idx+ii_x)=-1;
                        ub(st_idx+ii_x)=1;
                    end

                end                
                
                % fix y: 10,11
                if i_x>=10 && i_x<=11
                    if ii_x==2
                        lb(st_idx+ii_x)=-1;
                        ub(st_idx+ii_x)=1;
                    end

                end

                if i_x==2
                    if ii_x==2
                        lb(st_idx+ii_x)=-5;
                        ub(st_idx+ii_x)=20;
                    end
                end
            else
                % radius
                x0(st_idx+ii_x)=4;
                lb(st_idx+ii_x)=0;
                if i_x==1
                    lb(st_idx+ii_x)=2;
                    ub(st_idx+ii_x)=6;
                elseif i_x<=2
                    lb(st_idx+ii_x)=2;
                    ub(st_idx+ii_x)=6;
                elseif i_x<=5
                    ub(st_idx+ii_x)=15;
                elseif i_x>=7 && i_x<=9
                    ub(st_idx+ii_x)=25;
                elseif i_x==10
                    ub(st_idx+ii_x)=18;
                else
                    ub(st_idx+ii_x)=20;
                end
            end
                
        end
    end

    num_x0=size(x0,1);

    %----------------------------------fill cur_csg
    sketch_csg_cell=cell(num_x0,1);
    methods_cells=cell(num_x0,1);
    dsign_cell=cell(num_x0,1);    
    addtional_args=cell(num_x0,1);

    for i_x=1:num_cbs
        st_idx=(i_x-1)*num_c_in_cbs;
        for ii_x=1:num_c_in_cbs
            sketch_csg_cell{st_idx+ii_x}=skt_csg_list(i_x);
            dsign_cell{st_idx+ii_x}=-1;
            if ii_x==1
                methods_cells{st_idx+ii_x}='dx0';
                addtional_args{st_idx+ii_x}={};
            elseif ii_x==2
                methods_cells{st_idx+ii_x}='dy0';
                addtional_args{st_idx+ii_x}={}; 
            else
                methods_cells{st_idx+ii_x}='dRi';
                addtional_args{st_idx+ii_x}=ii_x-1; 
            end

        end
    end
    %---------------------------------fill deriv_cell
    deriv_cell=cell(num_x0,1);
    for i=1:num_x0
        sub_cell=cell(4,1);
        sub_cell{1}=sketch_csg_cell{i};
        sub_cell{2}=methods_cells{i};
        sub_cell{3}=dsign_cell{i};
        sub_cell{4}=addtional_args{i};
        deriv_cell{i}=sub_cell;
    end
    %---------------------------------end fill
    
    mp.handleMap('x0')=x0;
    mp.handleMap('lb')=lb;
    mp.handleMap('ub')=ub;
    mp.handleMap('modify_func')=@modify_c_by_x;
    mp.handleMap('constrast_list')= ...
        {};
    mp.handleMap('deriv_cell')=deriv_cell;
    mp.handleMap('target_rootH_list')=csg_list;
%     mp.handleMap('sketch_list')=skt_list;
    mp.handleMap('sketch_csg_list')=skt_csg_list;

    mp.handleMap('cbs_list')=cbs_list;
    modify_c_by_x(x0,mp)
end


function modify_c_by_x(x,mp)
    num_cbs=11;
    num_c_in_cbs=24;
    cbs_list=mp.handleMap('cbs_list');
    for i_x=1:num_cbs
        st_idx=(i_x-1)*num_c_in_cbs;
        for ii_x=1:num_c_in_cbs
            cur_cbs=cbs_list(i_x);
            if ii_x<=2
                cur_cbs.Oa(ii_x)=x(st_idx+ii_x);
            else
                cur_cbs.CVs(ii_x-1)=x(st_idx+ii_x);
            end

        end
    end    
end

function ret=cal_NeuF(X)
r=7.065;
h=7.62;
Or=[0 0];
fx_vec=[100 0 -100]/h*1e3*2/pi/r;
unit_fx=fx_vec([1,3])./norm(fx_vec([1,3]));

vec_X=X(:,[1,3])-Or;
unit_X=vec_X./vecnorm(vec_X,2,2);
the_cos_x=unit_X*unit_fx';

the_cos_x(the_cos_x<0)=0;

ret=unit_X.*(the_cos_x.*norm(fx_vec([1,3])) );

ret=[ret(:,1),zeros(size(ret,1),1),ret(:,2)];

end