function mp= Build_Piston28()
% Bulid mp obj （contains models and all arguments)
%% init
pace=3.33;
Pst=[-13 -20 -10]-pace*0.52;
Ped=[130 20 10];
Ped=FCM_dV.round_Ped(Pst,Ped,pace);
maxdeep=4;
order=2;
viewxy=[0 60];

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
%     SuiteSparse_paths2();
additon_path=pwd;

if count(py.sys.path,additon_path) == 0
    insert(py.sys.path,int32(0),additon_path);
end

%     additon_path=[additon_path,'/TopCSG'];
%     if count(py.sys.path,additon_path) == 0
%         insert(py.sys.path,int32(0),additon_path);
%     end




%% assemble c0
c=CSG("STL",[pwd,'/SolidWoks Files/Piston28/piston28.STL'],'init');
c1center=[25 0 0];
c2center=[55 0 0];

cir1=Circle(c1center,5);
cir2=Circle(c2center,5);

tLine1=TangentLine(cir1,cir2,'bot');
tLine2=TangentLine(cir1,cir2,'top');

tLine3=Line(tLine1.p2,tLine2.p1);
tLine4=Line(tLine2.p2,tLine1.p1);
tmp=tLine1*tLine3*tLine2*tLine4+cir1+cir2;

tLine1.addlistener('propChanged',@tLine4.syn_p21);
tLine1.addlistener('propChanged',@tLine3.syn_p12);
tLine2.addlistener('propChanged',@tLine3.syn_p21);
tLine2.addlistener('propChanged',@tLine4.syn_p12);
cir1.addlistener('propChanged',@tLine1.update_p);
cir1.addlistener('propChanged',@tLine2.update_p);
cir2.addlistener('propChanged',@tLine1.update_p);
cir2.addlistener('propChanged',@tLine2.update_p);

init_h=2;
%---------------------------------------top
the_skt=Sketch([0 0 init_h],Ra,5,tmp);
the_rc1=RoundCorner(3.5-init_h,the_skt,tLine1,false,true);
the_rc2=RoundCorner(3.5-init_h,the_skt,tLine2,false,true);
the_rcc1=RoundCorner_Circle(3.5-init_h,the_skt,cir1,false,true);
the_rcc2=RoundCorner_Circle(3.5-init_h,the_skt,cir2,false,true);

lines_skt=CSG();lines_skt.init(the_skt);
rc1_skt=CSG();rc1_skt.init(the_rc1,'Lrc1');
rc2_skt=CSG();rc2_skt.init(the_rc2,'Lrc2');
rcc1_skt=CSG();rcc1_skt.init(the_rcc1,'Crc1');
rcc2_skt=CSG();rcc2_skt.init(the_rcc2,'Crc2');
half_pla1=CSG('Plane3D',c1center+[0.1 0 0],[-1 0 0],'planeL');
half_pla2=CSG('Plane3D',c2center-[0.1 0 0],[1 0 0],'planeR');
%---------------------------------------down
the_skt_down=Sketch([0 0 -init_h],[pi 0 0],5,tmp);
the_rc1_down=RoundCorner(3.5-init_h,the_skt_down,tLine1,false,true);
the_rc2_down=RoundCorner(3.5-init_h,the_skt_down,tLine2,false,true);
the_rcc1_down=RoundCorner_Circle(3.5-init_h,the_skt_down,cir1,false,true);
the_rcc2_down=RoundCorner_Circle(3.5-init_h,the_skt_down,cir2,false,true);
lines_skt_down=CSG();lines_skt_down.init(the_skt_down);
rc1_skt_down=CSG();rc1_skt_down.init(the_rc1_down,'Lrc3');
rc2_skt_down=CSG();rc2_skt_down.init(the_rc2_down,'Lrc4');
rcc1_skt_down=CSG();rcc1_skt_down.init(the_rcc1_down,'Crc3');
rcc2_skt_down=CSG();rcc2_skt_down.init(the_rcc2_down,'Crc4');


c=c-lines_skt-lines_skt_down+(rc1_skt+rc1_skt_down+rc2_skt+rc2_skt_down+((rcc1_skt+rcc1_skt_down)*half_pla1)+((rcc2_skt+rcc2_skt_down)*half_pla2));

V0=1;
% time and calculate Volume
tic
fcmv=FCM_V(Pst,Ped,pace,maxdeep,order,@c.calSDF,@(X) ones(size(X,1),1));% |||1 dV
% V0=fcmv.calIntegral_V();%%
% fcmv=FCM_V(Pst,Ped,pace,maxdeep,order,@c.calSDF,@(X) ones(24,24,size(X,1)));% |||1 dV
% K0=sum(fcmv.cal_Ke(),3);
% K0(1)

% tmp_c=CSG('Sphere',Oa,5);
% BC=tmp_c.copy()*CSG('Plane3D',[0 0 0],[-1 0 0]);
% fcmdv=FCM_dV(Pst,Ped,pace,maxdeep,order,@c.calSDF,@(X) ones(size(X,1),1));% |||1 dV
% 
% S0=fcmdv.calIntegral_S(Dir)
% Sstd=4*pi*5*5;
% abs(S0-Sstd)/Sstd
% 
% fcmdv=FCM_V(Pst,Ped,pace,maxdeep,order,@tmp_c.calSDF,@(X) ones(size(X,1),1));% |||1 dV
% Vstd=4/3*pi*5^3;
% V00=fcmdv.calIntegral_V()
% abs(Vstd-V00)/Vstd
fprintf('V0: %f   ,Time: %f\n',V0,toc);

%% build Dirichlet/Neumann boundary
eepsD=1e-1*9;
eepsN=9e-2;
Dc1=Circle(Oa,9);
Dc2=Circle(Oa,9-eps);
Dcskt=(Dc1-Dc2);
skt=Sketch([0 0 -8],Ra,8*2,Dcskt);
Dirich_csg=CSG();
Dirich_csg.init(skt,'Dir');
Dirich_csg.judge_func=@(X) abs(vecnorm(X(:,1:2)-[0 0],2,2)-9)<eepsD & abs(X(:,3))<8+eepsD ;
% Dirich_csg.eeps=pace/8;
%------------------------------------------Neu
Nc1=Circle(Oa,15);
Nc2=Circle(Oa,15-eps);
% half_line=Line([0 -15 0],[0 15 0]);
Ncskt=(Nc1-Nc2);
skt=Sketch([110 0 -8],Ra,8*2,Ncskt);
Neuman_csg=CSG();
Neuman_csg.init(skt,'Neu');
Neuman_csg.judge_func=@(X) abs(vecnorm(X(:,1:2)-[110 0],2,2)-15)<eepsN & abs(X(:,3))<8+eepsN;

Neuman_csg.fdv=@cal_NeuF;
% Neuman_csg.eeps=Dirich_csg.eeps;

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=18*pi*18;
relaS=(S0-Sstd)/Sstd

S1=fcmdv.calIntegral_S(Neuman_csg)
S1std=2*pi*15*16;
relaS1=(S1-S1std)/S1std

%% show model
clf;
%---------------step1
c.set_display_args(Pst,Ped,pace/6,'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.delta=Dirich_csg.delta;
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

tmphandle=TmpHandle();
tmphandle.rcc1_skt=rcc1_skt;
tmphandle.rcc2_skt=rcc2_skt;
tmphandle.rc1_skt=rc1_skt;
tmphandle.rc2_skt=rc2_skt;
tmphandle.lines_skt=lines_skt;

tmphandle.rc1_skt_down=rc1_skt_down;
tmphandle.rc2_skt_down=rc2_skt_down;
tmphandle.rcc1_skt_down=rcc1_skt_down;
tmphandle.rcc2_skt_down=rcc2_skt_down;
tmphandle.lines_skt_down=lines_skt_down;

tmphandle.cir1=cir1;
tmphandle.cir2=cir2;
tmphandle.the_skt=the_skt;
tmphandle.the_skt_down=the_skt_down;

tmphandle.the_rc1=the_rc1;
tmphandle.the_rc2=the_rc2;
tmphandle.the_rcc1=the_rcc1;
tmphandle.the_rcc2=the_rcc2;

tmphandle.the_rc1_down=the_rc1_down;
tmphandle.the_rc2_down=the_rc2_down;
tmphandle.the_rcc1_down=the_rcc1_down;
tmphandle.the_rcc2_down=the_rcc2_down;
mp.tmphandle=tmphandle;

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

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

function ret=cal_NeuF(X)
r=15;
h=16;
Or=[110 0];
fx_vec=[-100 0 0]/h*1e3*2/pi/r;% 2/pi is coef 
fy_vec=[0 -200 0]/h*1e3*2/pi/r;
unit_fx=fx_vec(1:2)./norm(fx_vec(1:2));
unit_fy=fy_vec(1:2)./norm(fy_vec(1:2));

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

the_cos_x(the_cos_x<0)=0;
the_cos_y(the_cos_y<0)=0;

ret=unit_X.*(the_cos_x.*abs(fx_vec(1))+ the_cos_y.*abs(fy_vec(2)));

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

end