fem二维网格划分

Posted 2020-09-25 I know you

 主要考虑非结构网格

详细文档https://wenku.baidu.com/view/5185c6d680eb6294dd886c99.html?from=search

  

  

 

 

 

 

按照上面的思路分为如下几步:

测试主程序：

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)形成初始三角形

思路：1.扫描边界所有点，如果以该点及相邻两点形成的三角形具有最大的三角形最小内角，且没有其他边界点在这个三角形中，则进行连接

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