SciDog我先添加一段代码试试：
[cc lang=matlab]
classdef PicMesh < handle
properties
imgname
pixelsize
bw
m
n
end
properties(Dependent)
Ncell
NfaceX
NfaceY
Npoint
end
properties
ND = 2;
PFlag = 10;
FFlag = 13;
CFlag = 12;
NFlag = 2;
CommentFlag = 0;
end
properties
p
c
fx
fy
end
properties
Point = struct(‘index’,1,’type’,1,’data’,[]);
Face
Cell
end
methods
function obj = PicMesh(img)
bw = image_read(img);
bw = bw_filter(bw);
obj.bw = bwcrop(bw);
[m,n] = size(obj.bw);
obj.m = m;
obj.n = n;
obj.getCellIndex();
obj.getPoint();
obj.getFace();
end
end
methods
function obj = getCellIndex(obj)
[ci,cj] = find(obj.bw);
obj.c = [ci,cj];
P = cell(obj.Ncell,1);
Fx = cell(obj.Ncell,1);
Fy = cell(obj.Ncell,1);
for k = 1:obj.Ncell
[P{k},Fx{k},Fy{k}] = obj.genPointFace(ci(k),cj(k));
end
P = cell2mat(P);
P = unique(P,’row’);
Fx = cell2mat(Fx);
Fx = unique(Fx,’row’);
Fy = cell2mat(Fy);
Fy = unique(Fy,’row’);
obj.p = P;
obj.fx = Fx;
obj.fy = Fy;
end
function obj = getPoint(obj)
obj.Point.data = obj.p;
end
function obj = getFace(obj)
data = zeros(obj.NfaceX+obj.NfaceY,4);
boundarytype = zeros(obj.Ncell,1);
% data(:,1) = 2;
for k = 1:obj.NfaceX
[nt,nb,cl,cr] = obj.genFaceX(obj.fx(k,1),obj.fx(k,2));
N0 = obj.findInd(nt,obj.p);
N1 = obj.findInd(nb,obj.p);
CR = obj.findInd(cr,obj.c);
CL = obj.findInd(cl,obj.c);
data(k,:) = [N0 N1 CR CL];
if CR ==0
if obj.fx(k,2)==obj.n+1
boundarytype(k) = 7;
else
boundarytype(k) = 4;
end
end
if CL ==0
if obj.fx(k,2)==1
boundarytype(k) = 5;
else
boundarytype(k) = 4;
end
end
end
for k = 1:obj.NfaceY
[nl,nr,ct,cb] = obj.genFaceY(obj.fy(k,1),obj.fy(k,2));
N0 = obj.findInd(nl,obj.p);
N1 = obj.findInd(nr,obj.p);
CR = obj.findInd(ct,obj.c);
CL = obj.findInd(cb,obj.c);
data(obj.NfaceX+k,:) = [N0 N1 CR CL];
if CR ==0
if obj.fy(k,1) == 1
boundarytype(obj.NfaceX+k) = 6;
else
boundarytype(obj.NfaceX+k) = 4;
end
end
if CL ==0
if obj.fy(k,1) == obj.m+1
boundarytype(obj.NfaceX+k) = 8;
else
boundarytype(obj.NfaceX+k) = 4;
end
end
end
% reverse face
flag = data(:,3) == 0;
data(flag,:) = data(flag,[2 1 4 3]);
% 3 is interior
face(1).index = 3;
face(1).data = data(boundarytype == 0,:);
face(1).type = 2;
face(1).start = 1;
% 4 is wall but not boundary
face(2).index = 4;
face(2).data = data(boundarytype == 4,:);
face(2).type = 3;
face(2).start = face(1).start + numel(face(1).data(:,1));
% 5 is left boundary
face(3).index = 5;
face(3).data = data(boundarytype == 5,:);
face(3).type = 3;
face(3).start = face(2).start + numel(face(2).data(:,1));
% 6 is upper boundary
face(4).index = 6;
face(4).data = data(boundarytype == 6,:);
face(4).type = 3;
face(4).start = face(3).start + numel(face(3).data(:,1));
% 7 is right boundary
face(5).index = 7;
face(5).data = data(boundarytype == 7,:);
face(5).type = 3;
face(5).start = face(4).start + numel(face(4).data(:,1));
% 8 is bottom boundary
face(6).index = 8;
face(6).data = data(boundarytype == 8,:);
face(6).type = 3;
face(6).start = face(5).start + numel(face(5).data(:,1));
obj.Face = face;
end
function obj = getCell(obj)
obj.Cell.index = 2;
end
end
methods (Static)
function [P,Fx,Fy] = genPointFace(i,j)
P = [i,j;i+1,j;i,j+1;i+1,j+1];
Fx = [i,j;i,j+1];
Fy = [i,j;i+1,j];
end
function [nt,nb,cl,cr] = genFaceX(i,j)
nt = [i, j];
nb = [i+1,j];
cl = [i,j-1];
cr = [i, j];
end
function [nl,nr,ct,cb] = genFaceY(i,j)
nl = [i, j];
nr = [i,j+1];
ct = [i-1,j];
cb = [i, j];
end
function index = findInd(x,data)
flag = all(x == data,2);
if any(flag)
index = find(flag);
else
index = 0;
end
end
end
methods
function output(obj,filename)
fid = fopen(strcat(filename,’.msh’),’w’);
% first line with comment and company
fprintf(fid,'(0 “PIC2MSH to Fluent File, by SINOPEC:”)\r\n’);
fprintf(fid,’\r\n’);
% Dimension
fprintf(fid,'(0 “Dimension:”)\r\n’);
fprintf(fid,'(2 %d)\r\n’,obj.ND);
fprintf(fid,’\r\n’);
% Points
fprintf(fid,'(0 “Points:”)\r\n’);
fprintf(fid,'(10 (0 1 %X 1 2))\r\n’,obj.Npoint);
fprintf(fid,'(10 (1 1 %X 1 2)’,obj.Npoint);
fprintf(fid,'(\r\n’);
fprintf(fid,’%15.10e %15.10e \r\n’,obj.Point.data.’);
fprintf(fid,’))\r\n’);
fprintf(fid,’\r\n’);
% Faces
fprintf(fid,'(0 “Faces:”)\r\n’);
fprintf(fid,'(13 (0 1 %X 0))\r\n’,obj.NfaceX+obj.NfaceY);
for i = 1:numel(obj.Face)
fprintf(fid,'(13 (%d %x %x %d 2)’,…
obj.Face(i).index, obj.Face(i).start,…
obj.Face(i).start+numel(obj.Face(i).data(:,1))-1,…
obj.Face(i).type);
fprintf(fid,'(\r\n’);
fprintf(fid,’%x %x %x %x\r\n’,obj.Face(i).data.’);
fprintf(fid,’))\r\n’);
end
fprintf(fid,’\r\n’);
% Cell
fprintf(fid,'(0 “Cells:”)\r\n’);
fprintf(fid,'(12 (0 1 %X 0))\r\n’,obj.Ncell);
fprintf(fid,'(12 (2 1 %X 1 3))\r\n’,obj.Ncell);
fprintf(fid,’\r\n’);
% Zone
fprintf(fid,'(0 “Zones:”)\r\n’);
fprintf(fid,'(45 (2 fluid fluid)())\r\n’);
fprintf(fid,'(45 (3 interior interior-face)())\r\n’);
fprintf(fid,'(45 (4 wall wall.0)())\r\n’);
fprintf(fid,'(45 (5 wall wall.b)())\r\n’);
fprintf(fid,'(45 (6 wall wall.l)())\r\n’);
fprintf(fid,'(45 (7 wall wall.t)())\r\n’);
fprintf(fid,'(45 (8 wall wall.r)())\r\n’);
% finish file
fclose(fid);
end
end
methods
function Ncell = get.Ncell(obj)
Ncell = sum(obj.bw(:));
end
function NfaceX = get.NfaceX(obj)
NfaceX = numel(obj.fx(:,1));
end
function NfaceY = get.NfaceY(obj)
NfaceY = numel(obj.fy(:,1));
end
function Npoint = get.Npoint(obj)
Npoint = numel(obj.Point.data(:,1));
end
end
end
[/cc]
