You can not select more than 25 topics
Topics must start with a letter or number, can include dashes ('-') and can be up to 35 characters long.
222 lines
4.8 KiB
222 lines
4.8 KiB
function mp=Build_MMC()
|
|
% Bulid mp obj （contains models and all arguments)
|
|
%% init
|
|
pace=0.3;
|
|
Pst=[0 0 0]-pace*0.51;
|
|
Ped=[2 4 2];
|
|
Ped=FCM_dV.round_Ped(Pst,Ped,pace);
|
|
maxdeep=3;
|
|
order=1;
|
|
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
|
|
|
|
|
|
%% 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
|
|
% cube1=CSG("Cube","corner",[0 0 0],Ra,[12 12 1]);
|
|
% cube2=CSG("Cube","corner",[0 0 7],Ra,[12 12 1]);
|
|
% c=cube1+cube2;
|
|
|
|
cyObjs=[];
|
|
diffCoord=[
|
|
0 0 0;
|
|
1 0 0;
|
|
1 1 0;
|
|
0 1 0;
|
|
];
|
|
|
|
rotateAngles=[
|
|
-1 1 0;
|
|
-1 -1 0;
|
|
1 -1 0;
|
|
1 1 0;
|
|
].*(ones(4,1)*[pi/5 pi/4 0]);
|
|
lenCy=sqrt(3);
|
|
for i=0:1:1
|
|
for j=0:1:3
|
|
for k=0:1:1
|
|
for idx=1:4
|
|
cy=CSG("Cylinder",[i j k]+diffCoord(idx,:),rotateAngles(idx,:),0.1,lenCy);
|
|
cyObjs=[cyObjs cy];
|
|
end
|
|
end
|
|
end
|
|
end
|
|
firstFlg=1;
|
|
cys=[];
|
|
for cy=cyObjs
|
|
if firstFlg==1
|
|
firstFlg=0;
|
|
cys=cy;
|
|
else
|
|
cys=cys+cy;
|
|
end
|
|
end
|
|
limitedCube=CSG("Cube","corner",[0 0 0],Ra,[2 4 2]);
|
|
c=(limitedCube*cys);
|
|
|
|
|
|
|
|
%---------------------------assemble c
|
|
|
|
% show model
|
|
clf;
|
|
c.set_display_args([Pst(1) Pst(2) Pst(3)],Ped,pace/5,'viewxy',viewxy);
|
|
% c.display_csg();
|
|
% c.csg2stl('right');
|
|
|
|
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=5e-1;
|
|
|
|
% DC2d=Rect([0 0 0],[12 12 0],false);
|
|
% Dskt=Sketch([0 0 0],[0 0 0],1e-9,DC2d);
|
|
% Dirich_csg=CSG();Dirich_csg.init(Dskt);
|
|
Dirich_csg=CSG("Cube","corner",[0 0 0],Ra,[2 0 2]);
|
|
Dirich_csg.saveAs('Dir');
|
|
Dirich_csg.judge_func=@(X) abs(X(:,2))< eepsD*0.35;
|
|
%------------------------------------------Neu
|
|
|
|
% Ncir=Rect([0 0 0],[12 12 0],false);
|
|
% Nskt=Sketch([0 0 8],[0 0 0],1e-9,Ncir);
|
|
% Neuman_csg=CSG();
|
|
% Neuman_csg.init(Nskt,'Neu');
|
|
Neuman_csg=CSG("Cube","corner",[1-0.3 4 0-0.3],Ra,[1+0.3 4 0+0.3]);
|
|
Neuman_csg.judge_func=@(X) abs(X(:,2)-4)<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=12*12;
|
|
relaSD=(S0-Sstd)/Sstd
|
|
|
|
S1=fcmdv.calIntegral_S(Neuman_csg)
|
|
S1std=0.25
|
|
relaSN=(S1-S1std)/S1std
|
|
|
|
%% show model
|
|
clf;
|
|
c.set_display_args([Pst(1) 0 Pst(3)],Ped,pace/5,'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=0;
|
|
mp.fixedNeu=1;
|
|
mp.iscontrast=0;
|
|
mp.isrecord=1;
|
|
mp.frac=3;
|
|
mp.maxloop=30;
|
|
mp.interp_eeps=pace/(2^(maxdeep-1));
|
|
mp.interp_eeps=pi/2/200;
|
|
% mp.strongNBC=1;
|
|
% mp.Neu_st=[1 4 0]-[1 1 0]*0.01;
|
|
% mp.Neu_ed=[1 4 0]+[1 1 0]*0.01;
|
|
% mp.strongDBC=1;
|
|
%% save handle here
|
|
for i=1:length(cyObjs)
|
|
mp.handleMap(['cyObj' num2str(i)])=cyObjs(i);
|
|
end
|
|
mp.handleMap('lenObj')=length(cyObjs);
|
|
|
|
%% 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
|
|
x0=[];
|
|
lb=[];
|
|
ub=[];
|
|
for cyObj=cyObjs
|
|
% x0=[x0 cyObj.rootH.r cyObj.rootH.tx cyObj.rootH.ty];
|
|
x0=[x0 cyObj.rootH.x cyObj.rootH.y cyObj.rootH.tx cyObj.rootH.ty cyObj.rootH.r cyObj.rootH.h];
|
|
lb=[lb cyObj.rootH.x-1 cyObj.rootH.y-1 -pi/2 -pi/2 cyObj.rootH.r*0 0];
|
|
ub=[ub cyObj.rootH.x+1 cyObj.rootH.y+1 pi/2 pi/2 cyObj.rootH.r*1.5 10];
|
|
end
|
|
x0=x0';
|
|
lb=lb';
|
|
ub=ub';
|
|
%---------------------------------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')= ...
|
|
{};
|
|
modify_c_by_x(x0,mp);
|
|
|
|
end
|
|
|
|
|
|
function modify_c_by_x(x,mp)
|
|
Len=mp.handleMap('lenObj');
|
|
for i=1:Len
|
|
curObj=mp.handleMap(['cyObj' num2str(i)]);
|
|
curObj.rootH.x=x(6*i-5);
|
|
curObj.rootH.y=x(6*i-4);
|
|
curObj.rootH.tx=x(6*i-3);
|
|
curObj.rootH.ty=x(6*i-2);
|
|
curObj.rootH.r=x(6*i-1);
|
|
curObj.rootH.h=x(6*i);
|
|
if(curObj.rootH.r<0.2)
|
|
curObj.rootH.r=0;
|
|
end
|
|
if(curObj.rootH.h<0.2)
|
|
curObj.rootH.h=0;
|
|
end
|
|
end
|
|
end
|
|
|
|
function ret=cal_NeuF(X)
|
|
|
|
lenX=size(X,1);
|
|
% 1 N
|
|
ret=zeros(lenX,3);
|
|
ret(:,3)=1000;
|
|
|
|
end
|
|
|
|
|