function mp= Build_Fountain()
% Bulid mp obj （contains models and all arguments)
    %% init
    pace=2;
    Pst=[-30 -30 0]-pace*0.52;
    Ped=[30 30 15];
    Ped=FCM_dV.round_Ped(Pst,Ped,pace); 
    maxdeep=4;
    order=1;
    viewxy=[0 15];

    
    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


    %% 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
    rect_tot=Rect([4 0 0],[20 15 0]);

    obs_cvs=[(4:2:18)'];
    num_c_in_obs=length(obs_cvs);
    obs_cvs=[obs_cvs,obs_cvs*0+15];
    obs=OBS(obs_cvs,3,[4 35 0]);
%     obs.addlistener('propChanged',@obs.auto_convert_cvs);

%     tmp_ra=RotationZ([0 0 0],[0 0 0],obs);
%     pla=Plane3D([0 0 15],[0 0 -1]);
%     pla_csg=CSG();pla_csg.init(pla);
%     obs_csg=CSG();obs_csg.init(tmp_ra);
%     obs_csg=obs_csg*pla_csg;

    diff3=Diff2D(obs,-3);
    obs_margin=(diff3-obs)*Line([4 0 0],[4 -1 0])*Line([18 0 0],[18 1 0])+Rect([18 12 0],[20 15 0]);
    tmp_ra=RotationZ([0 0 0],[0 0 0],obs_margin);
    margin_csg=CSG();margin_csg.init(tmp_ra);

    
    num_c_in_cbs=10;
    cbs_cvs=zeros(num_c_in_cbs,1)+2;
    cbs_objs=[];%%
    for i=[4 12 20]
        for j=[0 10]
            tmp_cbs=CBS(cbs_cvs,3,[i j 0]);
            cbs_objs=[cbs_objs tmp_cbs];
        end
    end
    
    for i=[8 16]
        for j=[5]
            tmp_cbs=CBS(cbs_cvs,3,[i j 0]);
            cbs_objs=[cbs_objs tmp_cbs];
        end
    end
    
    %---------------------------deal each cbs_obj
    center_xy=zeros(0,2);
    cbs_csg_objs=[];%%
    for i_obj=1:length(cbs_objs)
        cur_obj=cbs_objs(i_obj);
        center_xy=cat(1,center_xy,cur_obj.Oa(1:2));

        tmp_ra=RotationZ([0 0 0],[0 0 0],cur_obj);
        tmp_CSG=CSG();tmp_CSG.init(tmp_ra);
        cbs_csg_objs=[cbs_csg_objs tmp_CSG];
    end
    %---------------------------
    
    c2d=(rect_tot-obs);
    
    %--------------------------------------------ears
    ear2d=Rect([-5 -25 0],[5 -19 0])+Circle([0 -25 0],5)-Circle([0 -25 0],2.5);
    ear_skt1=Sketch([0 0 15-3],[0 0 0],3,ear2d);      ear_c1=CSG();ear_c1.init(ear_skt1);
    ear_skt2=Sketch([0 0 15-3],[0 0 2*pi/3],3,ear2d); ear_c2=CSG();ear_c2.init(ear_skt2);
    ear_skt3=Sketch([0 0 15-3],[0 0 -2*pi/3],3,ear2d);ear_c3=CSG();ear_c3.init(ear_skt3);
    ear_c=ear_c1+ear_c2+ear_c3;
    
    
    [tx,ty]=meshgrid(-1:0.1:22);
    ret=c2d.calSDF([tx(:) ty(:) tx(:)*0]);
    figure
    contour(tx,ty,reshape(ret,size(tx)),ShowText="on");

[tx,ty,tz]=meshgrid(-30:0.3:30);
figure
the_x=[tx(:) ty(:) tz(:)];
ra=RotationZ([0 0 0],[0 0 0],c2d);
c=CSG();c.init(ra);

union_cbs_csg=cbs_csg_objs(1);
for i=2:length(cbs_csg_objs)
    union_cbs_csg=union_cbs_csg+cbs_csg_objs(i);
end

c=c-union_cbs_csg+margin_csg+ear_c;
ret=c.calSDF(the_x);
figure
isosurface(tx,ty,tz,reshape(ret,size(tx)),0);
axis equal

    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=1e-1;

DC2d=Circle([0 -25 0],2.5)-Circle([0 -25 0],2.5-eps);
DC_skt1=Sketch([0 0 15-3],[0 0 0],3,DC2d);      DC_c1=CSG();DC_c1.init(DC_skt1);
DC_skt2=Sketch([0 0 15-3],[0 0 2*pi/3],3,DC2d); DC_c2=CSG();DC_c2.init(DC_skt2);
DC_skt3=Sketch([0 0 15-3],[0 0 -2*pi/3],3,DC2d);DC_c3=CSG();DC_c3.init(DC_skt3);
Dirich_csg=DC_c1+DC_c2+DC_c3;
Dirich_csg.saveAs('Dir');
Dirich_csg.judge_func=@(X)  abs(X(:,1))<3 & abs(vecnorm(X(:,1:2)-[0 -25],2,2)-2.5)<eepsD | X(:,1)<-5 & abs(vecnorm(X(:,1:2)-[-21.65 12.50],2,2)-2.5)<eepsD |X(:,1)>5 & abs(vecnorm(X(:,1:2)-[21.65 12.50],2,2)-2.5)<eepsD;
%------------------------------------------Neu

obs_out=Diff2D(obs,0+eps);
Nc2d=(obs_out-obs)*Line([4 0 0],[4 -1 0])*Line([18 0 0],[18 1 0]);
NcRa=RotationZ([0 0 0],[0 0 0],Nc2d);
Neuman_csg=CSG();
Neuman_csg.init(NcRa,'Neu');
Neuman_csg.judge_func=[];
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=2*pi*2.5*3 *3;
relaSD=(S0-Sstd)/Sstd

S1=fcmdv.calIntegral_S(Neuman_csg)
S1std=pi*(18^2-4^2);
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=0;
    mp.iscontrast=0;
    mp.isrecord=0;
    mp.frac=0.6;
    mp.maxloop=60;
    mp.interp_eeps=pace/(2^(maxdeep-1));
    mp.Neu_norm_flg=1;
    %% save handle here


    mp.handleMap('obs')=obs;
    mp.handleMap('cbs_list')=cbs_objs;
    
    %% 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
    % 1-7:obs's y axis;  8-23:cbs*8's xy axis; 24-87:cbs*8's r
    x0=[zeros(1,num_c_in_obs-1)+15, center_xy(:)',                 zeros(1,(num_c_in_cbs-2)*8)+2 ]';
    lb=[zeros(1,num_c_in_obs-1)+5,  zeros(1,8)-5,zeros(1,8)-5,    zeros(1,(num_c_in_cbs-2)*8)+0 ]';
    ub=[zeros(1,num_c_in_obs-1)+15, zeros(1,8)+15,zeros(1,8)+20   , zeros(1,(num_c_in_cbs-2)*8)+15]';
    

    num_x0=size(x0,1);
    %--------------------------------fill target_rootH_cell
    target_rootH_cell=cell(num_x0,1);
    for i=1:7
        target_rootH_cell{i}=Neuman_csg;
    end
    cnt=0;
    for i=8:15
        cnt=cnt+1;
        target_rootH_cell{i}=cbs_csg_objs(cnt);
    end
    cnt=0;
    for i=16:23
        cnt=cnt+1;
        target_rootH_cell{i}=cbs_csg_objs(cnt);
    end

    for i=24:87
        cur_i=floor((i-24)/(num_c_in_cbs-2))+1;
        target_rootH_cell{i}=cbs_csg_objs(cur_i);
    end
    %--------------------------------fill changed_rootH_cell
    changed_rootH_cell=cell(num_x0,1);
    for i=1:7
        changed_rootH_cell{i}=obs;
    end
    cnt=0;
    for i=8:15
        cnt=cnt+1;
        changed_rootH_cell{i}=cbs_objs(cnt);
    end
    cnt=0;
    for i=16:23
        cnt=cnt+1;
        changed_rootH_cell{i}=cbs_objs(cnt);
    end

    for i=24:87
        cur_i=floor((i-24)/(num_c_in_cbs-2))+1;
        changed_rootH_cell{i}=cbs_objs(cur_i);
    end
    %--------------------------------fill dsign_cell
    dsign_cell=cell(num_x0,1);
    for i=1:87
        dsign_cell{i}=-1;
    end
    %--------------------------------fill prop_cell
    props_cell=cell(num_x0,1);
    for i=1:7
        props_cell{i}='CVs0';
    end
    for i=8:23
        props_cell{i}='Oa';
    end

    for i=24:87
        props_cell{i}='CVs';
    end
    %--------------------------------fill prop_idx_cell
    props_idx_cell=cell(num_x0,1);
    for i=1:7
        props_idx_cell{i}=8+i;
    end
    for i=8:15
        props_idx_cell{i}=1;
    end
    for i=16:23
        props_idx_cell{i}=2;
    end
    for i=24:87
        cur_i=mod(i-24,num_c_in_cbs-2)+1+1;
        props_idx_cell{i}=cur_i;
    end    
    %---------------------------------fill deriv_cell
    deriv_cell=cell(num_x0,1);
    for i=1:num_x0
        sub_cell=cell(5,1);
        sub_cell{1}=target_rootH_cell{i};
        sub_cell{2}=changed_rootH_cell{i};
        sub_cell{3}=dsign_cell{i};
        sub_cell{4}=props_cell{i};
        sub_cell{5}=props_idx_cell{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')=[cbs_csg_objs Neuman_csg];
    modify_c_by_x(x0,mp)
end


function modify_c_by_x(x,mp)
    num_c_in_obs=8;
    num_c_in_cbs=10;
    obs=mp.handleMap('obs');
    obs.CVs0(1:(end-1),2)=x(1:num_c_in_obs-1);
    obs.CVs=obs.convert_cvs(obs.CVs0,obs.Oa);

    cbs_list=mp.handleMap('cbs_list');
    for i_obj=1:length(cbs_list)
        cur_cbs=cbs_list(i_obj);
         cur_cbs.Oa(1)=x(num_c_in_obs-1   +i_obj);
         cur_cbs.Oa(2)=x(num_c_in_obs-1+8 +i_obj);
         cur_cbs.CVs(2:end-1)=x(num_c_in_obs-1+16  +(i_obj-1)*(num_c_in_cbs-2)+(1:num_c_in_cbs-2));
    end
end

function ret=cal_NeuF(X,norm_X)
%     ret=norm_X*1e3;
    ret=zeros(size(X));
    ret(:,3)=-1e3;
end