# fem二维网格划分

主要考虑非结构网格

clc;
%clear all ;
addpath(genpath(pwd));
pos=[0,5;0,0;2,0;2,2; ...
3,2;3,0;8,0;8,5];
boundary=[1 2 3 4 5 6 7 8 1];
length_max=0.5;
[pos,nie]=makeGrid(pos,boundary,length_max);


划分网格函数主程序

function [pos,nie]=makeGrid(pos,boundary,length_max)
%this function is to generate Grids for a single object
%boundary is a matrix describe the point-index of the boundary,
%described in a counter-clock direction.

%find the special points
[pos,boundary]=generate_Boundary_Points(pos,boundary,length_max);
% drawGrid(pos,nie);
[pos,nie]=initial_triangle(pos,boundary);
[pos,nie]=lop_all(pos,nie);
drawGrid(pos,nie);
[pos,nie]=fill_points(pos,nie,boundary,[0.2,2]);
end


（1）形成边界

本次程序采用均匀分割，设定一个最大的边界段值，如果代分割的边界比它大，则取中点

function [pos,boundary]=generate_Boundary_Points(pos,boundary,length_max)
iP=length(pos(:,1));
i=1;
while i<length(boundary)
if get_distance(pos,boundary(i),boundary(i+1))>length_max
pos=[pos;(pos(boundary(i),1)+pos(boundary(i+1),1))/2,(pos(boundary(i),2)+pos(boundary(i+1),2))/2];
iP=iP+1;%new Point
boundary=[boundary(1:i), iP,boundary(i+1:end)];
else
i=i+1;
end
end
end


(2)形成初始三角形

2.删除该边界点，更新边界；

3.重复，直到所有点处理完成。

function [pos,nie]=initial_triangle(pos,boundary)
nie=[];
boundary=[boundary(end-1),boundary];

while length(boundary)>4
i=2;
max_min_angle=get_min_angle(pos(boundary(i-1:i+1),:));
%get the i of the max of the min angle
for ii=3:length(boundary)-1
max_min_angle_temp=get_min_angle(pos(boundary(ii-1:ii+1),:));
if max_min_angle_temp>max_min_angle && ~is_any_inside_triangle(pos,boundary(ii-1:ii+1))
max_min_angle=max_min_angle_temp;
i=ii;
end
end
angle=get_point_angle(pos(boundary(i-1:i+1),:));

nie=[nie;boundary(i-1:i+1)];
drawGrid(pos,nie);
boundary(i)=[];
boundary(end)=boundary(2);
boundary(1)=boundary(end-1);
end


（3）　填充内部点函数

1.三角形最小内边，在所有三角形最大

2.试探新点位于该三角形的外心，且至少有一个相邻三角形外接圆包含该新点

3.填充一个新点后进行 局部优化（lop）和平滑操作（smooth）。

function [pos,nie]=fill_points(pos,nie,boundary,range)
flag=1;

while flag
flag=0;
max_edge=0;
%get the i of the max of the min angle
edge=[];
tri_delete=[];
for ii=1:length(nie)
[x_temp,y_temp,r]=get_outer_circle(pos(nie(ii,:),:));
[edge_temp,tri_delete_temp]=get_ploygon(pos,nie,ii,[x_temp,y_temp],0);
edges=get_edges(pos(nie(ii,:),:));
if  max_edge<min(edges)&& min(edges)>range(1) && length(edge_temp)>4
max_edge=min(edges);
edge=edge_temp;
tri_delete=tri_delete_temp;
x=x_temp;y=y_temp; %center point
ie=ii;
flag=1;
end

end
if mod(length(pos(:,1)),50)==0
flag=1;
end

%min_angle=get_min_angle(pos(nie(ie,:),:));

%[x,y]=get_new_point(pos(nie(ie,:),:));
[x,y,r]=get_outer_circle(pos(nie(ie,:),:));
% drawGrid(pos,nie);hold on
%   plot(x,y,'--rs');hold off

if length(edge)>4 && is_inside_polygon(pos,[x,y],edge)
tri_delete=[tri_delete,ie];
nie(tri_delete,:)=[]; %delete current element
pos=[pos;x,y];% add new point
ip_new=length(pos(:,1));
for i_edge=1:length(edge)-1;%add new elements
nie=[nie;ip_new,edge(i_edge),edge(i_edge+1)];
%  drawGrid(pos,nie);
end
end
[pos,nie]=lop_all(pos,nie);drawGrid(pos,nie);
pos=smooth_all(pos,nie,boundary);drawGrid(pos,nie);

end


posted @ 2017-06-26 23:04  I know you  阅读(691)  评论(0编辑  收藏