精华内容
下载资源
问答
  • 划分naca0012网格,其中interfunction为翼型函数
  • Matlab网格划分

    2012-10-29 21:36:21
    Matlab进行有限元的建模,并对有限元模型进行网格划分
  • matlab网格划分

    2013-06-08 19:42:53
    matlab网格划分,互相学习,但愿有帮助
  • MATLAB网格划分

    2010-07-24 17:35:30
    MATLAB语言划分网格,计算流体的基础
  • Matlab网格划分程序Distmesh讲解(一)

    万次阅读 2015-09-28 10:56:07
    Distmesh 是一个matlab语言写的网格划分软件。 源文件可以从上面的网址获取。 这里按行讲解各个算例。 p01_demo: 概算例是一个单位圆(半径为1)的网格划分,划分后的网格为: 以下逐行讲解该算例: ...
    http://persson.berkeley.edu/distmesh/
    

    Distmesh 是一个matlab语言写的网格划分软件。
    源文件可以从上面的网址获取。
    这里按行讲解各个算例。
    p01_demo:
    概算例是一个单位圆(半径为1)的网格划分,划分后的网格为:
    Matlab网格划分程序Distmesh讲解(一)
    以下逐行讲解该算例:
    function p01_demo ( iteration_max, h )
     Parameters:
     Input, integer ITERATION_MAX, the maximum number of iterations that DISTMESH
         should take.   (The program might take fewer iterations if it detects convergence.)
     Input, real h, the mesh spacing parameter.
    % 这里有两个输入参数,一个是ITERATION_MAX,迭代的最大次数。
    % 另一个是h, 网格划分的大小。 0<h<1
    % 默认参数值为: ITERATION=200     h=0.1


    [ p, t ]=distmesh_2d( fd, fh, h0, box, iteration_max, fixed );
    函数需要至少六个参数。
    d = fd ( p ),     p=[x y]
    fd   给定任一点到边界的距离函数,本例中定义为: d =   sqrt(x^2+y^2)-1;
    fh, scaled edge length function h(x,y). 也就是网格大小的函数。

    h0 也就是h, 网格的大小
    real BOX(2,2), the bounding box [xmin,ymin; xmax,ymax].
    最大外围矩形范围。 本例中为[0,0;1,1]
    ITERATION_MAX, the maximum number of iterations.
    real PFIX(NFIX,2), the fixed node positions. 网格中需要固定的点坐标,也就是一定需要出现在网格中的点。

    输出参数:
    real P(N,2), the node positions.   网格点的x,y坐标
    integer T(NT,3), the triangle indices. 输出网格任一一个三角形的三个顶点。

    第一步:
       [ x, y ] = meshgrid ( box(1,1) : h0                     : box(2,1), ...
                                                   box(1,2) : h0*sqrt(3)/2 : box(2,2) );
    根据h0,网格的大小,先把能涵盖欲划分区域的最大矩形划分为结构网格。
    然后把偶数行的点整体向右平移半格,
      x(2:2:end,:) = x(2:2:end,:) + h0 / 2;
    效果如下:
    Matlab网格划分程序Distmesh讲解(一)

    Matlab网格划分程序Distmesh讲解(一)
    第二步:
    根据fd的函数定义,移除边界外的点。
      p = p( feval_r( fd, p, varargin{:} ) <= geps, : );
    varagin为fd,fh的附加参数,这里为空。
    geps = 0.001 * h0;
    也就是保留了到边界的距离以外0.001 * h0以内的点。
    Matlab网格划分程序Distmesh讲解(一) 
    根据网格密度函数fh,每个点上产生一个0-1随机数,判断是否小于r0/max(r0)
    大于的话,改点被删除。
      p = [ pfix; p(rand(size(p,1),1) < r0 ./ max ( r0 ),: ) ];
       [ nfix, dummy ] = size ( pfix );

    当指定了某些点要保留的时候,把保留的点加入,删除重复的点。
     Especially when the user has included fixed points, we may have a few
     duplicates.   Get rid of any you spot.
    %
       p = unique ( p, 'rows' );
       N = size ( p, 1 );
    这个时候产生的网格如下:
    Matlab网格划分程序Distmesh讲解(一)

    第三步:迭代
       pold = inf; %第一次迭代前设置旧的点的坐标为无穷
    while ( iteration < iteration_max ) 
           iteration = iteration + 1;

    %先判断上次移动后的点和旧的点之间的移动距离,如果小于某个阀值,停止迭代
           if ( ttol < max ( sqrt ( sum ( ( p - pold ).^2, 2 ) ) / h0 ) )
     
               pold = p;     %如果还可以移动,保存当前的节点
               t = delaunayn ( p );   %利用delauny算法,生成三角形网格
               triangulation_count = triangulation_count + 1;
               pmid = ( p(t(:,1),:) + p(t(:,2),:) + p(t(:,3),:) ) / 3;   %计算三角形的重心。
               t = t( feval_r( fd, pmid, varargin{:} ) <= -geps, : ); % 移除重心在边界外部的三角形

    %        4. Describe each bar by a unique pair of nodes.
    %
    % 生成网格的边的集合,也就是相邻点之间连接的线段  
         bars = [ t(:,[1,2]); t(:,[1,3]); t(:,[2,3]) ];
               bars = unique ( sort ( bars, 2 ), 'rows' );

       end


      %
     6. Move mesh points based on bar lengths L and forces F
    %
     Make a list of the bar vectors and lengths.
     Set L0 to the desired lengths, F to the scalar bar forces,
     and FVEC to the x, y components of the bar forces.
    %
     At the fixed positions, reset the force to 0.
    %
           barvec = p(bars(:,1),:) - p(bars(:,2),:);     % 生成bar的矢量
           L = sqrt ( sum ( barvec.^2, 2 ) );   %计算bar的长度

       %根据每个bar的中点坐标,计算需要的三角形边的边长(这个在fh函数里控制)
           hbars = feval_r( fh, (p(bars(:,1),:)+p(bars(:,2),:))/2, varargin{:} );

    % 计算 需要的bar的长度,已经乘上了两个scale参数 Fscale, sqrt ( sum(L.^2) / sum(hbars.^2) );
    % 具体可参考他们的paper
           L0 = hbars * Fscale * sqrt ( sum(L.^2) / sum(hbars.^2) );

    % 计算每个bar上力
           F = max ( L0 - L, 0 );

    �r上力的分量,x,y方向                                                          
           Fvec = F ./ L * [1,1] .* barvec; 

    % 计算Ftot, 每个节点上力的残量                                  
           Ftot = full ( sparse(bars(:,[1,1,2,2]),ones(size(F))*[1,2,1,2],[Fvec,-Fvec],N,2) );

    %对于固定点,力的残量为零
           Ftot(1:size(pfix,1),:) = 0;

    % 根据每个节点上的受力,移动该点                                              
           p = p + deltat * Ftot; 

     7. Bring outside points back to the boundary
    %
     Use the numerical gradient of FD to project points back to the boundary.
    %

       
           d = feval_r( fd, p, varargin{:} ); %计算点到边界距离
           ix = d > 0;
         
           %计算移动梯度,相对边界
           dgradx = ( feval_r(fd,[p(ix,1)+deps,p(ix,2)],varargin{:}) - d(ix) ) / deps;
           dgrady = ( feval_r(fd,[p(ix,1),p(ix,2)+deps],varargin{:}) - d(ix) ) / deps;

         %将这些移动到边界外的投射回边界上
           p(ix,:) = p(ix,:) - [ d(ix) .* dgradx, d(ix) .* dgrady ];
    %
     I needed the following modification to force the fixed points to stay.
     Otherwise, they can drift outside the region and be lost.
     JVB, 21 August 2008.
    %
           p(1:nfix,1:2) = pfix;

           N = size ( p, 1 );
    %
     8. Termination criterion: All interior nodes move less than dptol (scaled)
    %
           if ( max ( sqrt ( sum ( deltat * Ftot ( d < -geps,:).^2, 2 ) ) / h0 ) < dptol )
               break;
           end

       end



    展开全文
  • matlab网格划分程序与matlab有限元的结合

    万次阅读 多人点赞 2014-02-11 11:10:42
    1. distmesh是一个较好的网格划分程序,具体可以参考:http://persson.berkeley.edu/distmesh/   2.matlab有限元可以参考徐荣桥的书 ...网格划分:  >>fh=@(p) 0.05+0.3*dellipse(p,[0.5,0.2]); >> fd=@(p) d

     1.

    distmesh是一个较好的网格划分程序,具体可以参考:http://persson.berkeley.edu/distmesh/

     

    2.matlab有限元可以参考徐荣桥的书

    3.这里本人打算画一个园内包含一个椭圆的模型:

     

    具体程序如下:

    a.

    网格划分:

     >>fh=@(p) 0.05+0.3*dellipse(p,[0.5,0.2]);
    >> fd=@(p) ddiff(dcircle(p,0,0,1),dellipse(p,[0.5,0.2]));
    >> [p,t]=distmesh2d(fd,fh,0.05,[-1,-1;1,1],[-1,-1;-1,1;1,-1;1,1]);

     

    b.在distmesh2d.m最后插入语句,导出数据格式(节点坐标,单元材料属性,边界约束条件)

    fid = fopen('exam4_2.dat3', 'w') ;
    [ilength,jlength]=size( p );
          fprintf(fid,'%d\n',ilength);
          for i=1:1:ilength
              fprintf(fid,'%d\t%f\t%f\n',i,p(i,1),p(i,2));
          end
          [ilength,jlength]=size( t );
          fprintf(fid,'%d\n',ilength);
          for i=1:1:ilength
              fprintf(fid,'%d\t%d\t%d\t%d\t%d\n',i,t(i,1),t(i,2),t(i,3),1);
          end
     fclose(fid);
    文件'exam4_2.dat3内容如下:

    5105(节点个数)
    1 -0.707106 -0.707106 (每个节点坐标)
    2 -0.707106 0.707107
    ...

    9869(三角单元个数)
    1 5006 4951 4934 1  (标号,i,j,k,材料属性的编号)
    2 197 147 100 1
    3 413 2 366 1
    4 413 2 323 1

    ...

    5(5种材料属性)
    1   2.00e9    0.25  100.0  22.0e3(编号,杨氏模量,泊松比,厚度[平面应力填1],密度)
    2   2.60e9    0.20  1.0  23.0e3
    3   2.85e10   0.20  1.0  25.0e3
    4   1.85e10   0.20  1.0  23.0e3
    5   2.85e10   0.20  1.0  22.0e3
    2(约束个数)
      1 2  0.0 (节点,方向【2,为y方向】,0.0(0为固定))
      1 1  0.0

      4 2  0.0 (节点,方向【2,为y方向】,0.0(0为固定))
      41  0.0

    c。调用有限元程序exam_2.m(见徐荣桥书)

    >>exam4_2 exam4_2.dat3

    d.后处理

    >> exam4_2_post

    最大主应力云图如下:

    最大剪应力云图如下:

     

    最小主应力云图如下:

     

    exam_2。m源程序如下:

    function exam4_2( file_in )
    % 本程序为第四章的第二个算例,采用三角形单元计算隧道在自重作用下的变形和应力
    %      exam4_2( filename )
    %  输入参数:
    %      file_in  ---------- 有限元模型文件

    % 定义全局变量
    %      gNode ------------- 节点坐标
    %      gElement ---------- 单元定义
    %      gMaterial --------- 材料性质
    %      gBC1 -------------- 第一类约束条件
    %      gK ---------------- 整体刚度矩阵
    %      gDelta ------------ 整体节点坐标
    %      gNodeStress ------- 节点应力
    %      gElementStress ---- 单元应力
        global gNode gElement gMaterial gBC1 gK gDelta gNodeStress gElementStress gNF

        if nargin < 1
            file_in = 'exam4_2.dat' ;
        end
       
        % 检查文件是否存在
        if exist( file_in ) == 0
            disp( sprintf( '错误:文件 %s 不存在', file_in ) )
            disp( sprintf( '程序终止' ) )
            return ;
        end
       
        % 根据输入文件名称生成输出文件名称
        [path_str,name_str,ext_str] = fileparts( file_in ) ;
        ext_str_out = '.mat' ;
        file_out = fullfile( path_str, [name_str, ext_str_out] ) ;
       
        % 检查输出文件是否存在
        if exist( file_out ) ~= 0
            answer = input( sprintf( '文件 %s 已经存在,是否覆盖? ( Yes / [No] ):  ', file_out ), 's' ) ;
            if length( answer ) == 0
                answer = 'no' ;
            end
           
            answer = lower( answer ) ;
            if answer ~= 'y' | answer ~= 'yes'
                disp( sprintf( '请使用另外的文件名,或备份已有的文件' ) ) ;
                disp( sprintf( '程序终止' ) );
                return ;
            end
        end

        % 建立有限元模型并求解,保存结果
        FemModel( file_in ) ;          % 定义有限元模型
        SolveModel ;                   % 求解有限元模型
        SaveResults( file_out ) ;      % 保存计算结果
       
        % 计算结束
        disp( sprintf( '计算正常结束,结果保存在文件 %s 中', file_out ) ) ;
        disp( sprintf( '可以使用后处理程序 exam4_2_post.m 显示计算结果' ) ) ;
    return ;

    function FemModel(filename)
    %  定义有限元模型
    %  输入参数:
    %      filename --- 有限元模型文件
    %  返回值:
    %      无
    %  说明:
    %      该函数定义平面问题的有限元模型数据:
    %        gNode ------- 节点定义
    %        gElement ---- 单元定义
    %        gMaterial --- 材料定义,包括弹性模量,梁的截面积和梁的抗弯惯性矩
    %        gBC1 -------- 约束条件

        global gNode gElement gMaterial gBC1 gNF
       
        % 打开文件
        fid = fopen( filename, 'r' ) ;
       
        % 读取节点坐标
        node_number = fscanf( fid, '%d', 1 ) ;
        gNode = zeros( node_number, 2 ) ;
        for i=1:node_number
            dummy = fscanf( fid, '%d', 1 ) ;
            gNode( i, : ) = fscanf( fid, '%f', [1, 2] ) ;
        end
       
        % 读取单元定义
        element_number = fscanf( fid, '%d', 1 ) ;
        gElement = zeros( element_number, 4 ) ;
        for i=1:element_number
            dummy = fscanf( fid, '%d', 1 ) ;
            gElement( i, : ) = fscanf( fid, '%d', [1, 4] ) ;
        end
       
        % 读取材料信息
        material_number = fscanf( fid, '%d', 1 ) ;
        gMaterial = zeros( material_number, 4 ) ;
        for i=1:material_number
            dummy = fscanf( fid, '%d', 1 ) ;
            gMaterial( i, : ) = fscanf( fid, '%f', [1,4] ) ;
        end
       
        % 读取边界条件
        bc1_number = fscanf( fid, '%d', 1 ) ;
        gBC1 = zeros( bc1_number, 3 ) ;
        for i=1:bc1_number
            gBC1( i, 1 ) = fscanf( fid, '%d', 1 ) ;
            gBC1( i, 2 ) = fscanf( fid, '%d', 1 ) ;
            gBC1( i, 3 ) = fscanf( fid, '%f', 1 ) ;
        end
       
        % 关闭文件
        fclose( fid ) ;
         % 集中力
        %     节点号   自由度号   集中力值
        gNF = [  1,       1,         -800e3;
           1,       2,         -800e3;
           4,       1,         -800e3;
           4,       2,         -800e3;
               ] ;
     
       
    return

    function SolveModel
    %  求解有限元模型
    %  输入参数:
    %     无
    %  返回值:
    %     无
    %  说明:
    %      该函数求解有限元模型,过程如下
    %        1. 计算单元刚度矩阵,集成整体刚度矩阵
    %        2. 计算单元的等效节点力,集成整体节点力向量
    %        3. 处理约束条件,修改整体刚度矩阵和节点力向量
    %        4. 求解方程组,得到整体节点位移向量
    %        5. 计算单元应力和节点应力

        global gNode gElement gMaterial gBC1 gK gDelta gNodeStress gElementStress gNF

        % step1. 定义整体刚度矩阵和节点力向量
        [node_number,dummy] = size( gNode ) ;
        gK = sparse( node_number * 2, node_number * 2 ) ;
        f = sparse( node_number * 2, 1 ) ;

        % step2. 计算单元刚度矩阵,并集成到整体刚度矩阵中
        [element_number,dummy] = size( gElement ) ;
        for ie=1:1:element_number
            disp( sprintf(  '计算刚度矩阵,当前单元: %d', ie  ) ) ;
            k = StiffnessMatrix( ie ) ;
            AssembleStiffnessMatrix( ie, k ) ;
        end

     % step3. 把集中力直接集成到整体节点力向量中
        [nf_number, dummy] = size( gNF ) ;
       for inf=1:1:nf_number
           n = gNF( inf, 1 ) ;
     
         d = gNF( inf, 2 ) ;
           f( (n-1)*2 + d ) = gNF( inf, 3 ) ;
      end
     
        % step3. 计算自重产生的等效节点力
    %     for ie=1:1:element_number
    %         disp( sprintf(  '计算自重的等效节点力,当前单元: %d', ie  ) ) ;
    %         egf = EquivalentGravityForce( ie ) ;
    %         i = gElement( ie, 1 ) ;
    %         j = gElement( ie, 2 ) ;
    %         m = gElement( ie, 3 ) ;
    %         f( (i-1)*2+1 : (i-1)*2+2 ) = f( (i-1)*2+1 : (i-1)*2+2 ) + egf( 1:2 ) ;
    %         f( (j-1)*2+1 : (j-1)*2+2 ) = f( (j-1)*2+1 : (j-1)*2+2 ) + egf( 3:4 ) ;
    %         f( (m-1)*2+1 : (m-1)*2+2 ) = f( (m-1)*2+1 : (m-1)*2+2 ) + egf( 5:6 ) ;
    %     end
     
        % step4. 处理约束条件,修改刚度矩阵和节点力向量。采用乘大数法
        [bc_number,dummy] = size( gBC1 ) ;
        for ibc=1:1:bc_number
            n = gBC1(ibc, 1 ) ;
            d = gBC1(ibc, 2 ) ;
            m = (n-1)*2 + d ;
            f(m) = gBC1(ibc, 3)* gK(m,m) * 1e15 ;
            gK(m,m) = gK(m,m) * 1e15 ;
        end

        % step 5. 求解方程组,得到节点位移向量
        gDelta = gK \ f ;
       
        % step 6. 计算单元应力
        gElementStress = zeros( element_number, 6 ) ;
        for ie=1:element_number
            disp( sprintf(  '计算单元应力,当前单元: %d', ie  ) ) ;
            es = ElementStress( ie ) ;
            gElementStress( ie, : ) = es ;
        end
       
        % step 7. 计算节点应力(采用绕节点加权平均)
        gNodeStress = zeros( node_number, 6 ) ;      
        for i=1:node_number                        
            disp( sprintf(  '计算节点应力,当前结点: %d', i  ) ) ;
            S = zeros( 1, 3 ) ;                        
            A = 0 ;
            for ie=1:1:element_number
                for k=1:1:3
                    if i == gElement( ie, k )
                        area= ElementArea( ie ) ;
                        S = S + gElementStress(ie,1:3 ) * area ;
                        A = A + area ;
                        break ;
                    end
                end
            end
            gNodeStress(i,1:3) = S / A ;
            gNodeStress(i,6) = 0.5*sqrt( (gNodeStress(i,1)-gNodeStress(i,2))^2 + 4*gNodeStress(i,3)^2 ) ;
            gNodeStress(i,4) = 0.5*(gNodeStress(i,1)+gNodeStress(i,2)) + gNodeStress(i,6) ;
            gNodeStress(i,5) = 0.5*(gNodeStress(i,1)+gNodeStress(i,2)) - gNodeStress(i,6) ;
        end
    return

    function k = StiffnessMatrix( ie )
    %  计算单元刚度矩阵
    %  输入参数:
    %     ie ----  单元号
    %  返回值:
    %     k  ----  单元刚度矩阵

        global gNode gElement gMaterial
        k = zeros( 6, 6 ) ;
        E  = gMaterial( gElement(ie, 4), 1 ) ;
        mu = gMaterial( gElement(ie, 4), 2 ) ;
        h  = gMaterial( gElement(ie, 4), 3 ) ;
        xi = gNode( gElement( ie, 1 ), 1 ) ;
        yi = gNode( gElement( ie, 1 ), 2 ) ;
        xj = gNode( gElement( ie, 2 ), 1 ) ;
        yj = gNode( gElement( ie, 2 ), 2 ) ;
        xm = gNode( gElement( ie, 3 ), 1 ) ;
        ym = gNode( gElement( ie, 3 ), 2 ) ;
        ai = xj*ym - xm*yj ;
        aj = xm*yi - xi*ym ;
        am = xi*yj - xj*yi ;
        bi = yj - ym ;
        bj = ym - yi ;
        bm = yi - yj ;
        ci = -(xj-xm) ;
        cj = -(xm-xi) ;
        cm = -(xi-xj) ;
        area = (ai+aj+am)/2 ;
        B = [bi  0 bj  0 bm  0
              0 ci  0 cj  0 cm
             ci bi cj bj cm bm] ;
        B = B/2/area ;
        D = [ 1-mu    mu      0
               mu    1-mu     0
                0      0   (1-2*mu)/2] ;
        D = D*E/(1-2*mu)/(1+mu) ;
        k = transpose(B)*D*B*h*abs(area) ;   
    return

    function B = MatrixB( ie )
    %  计算单元的应变矩阵B
    %  输入参数:
    %     ie ----  单元号
    %  返回值:
    %     B  ----  单元应变矩阵
        global gNode gElement
        xi = gNode( gElement( ie, 1 ), 1 ) ;
        yi = gNode( gElement( ie, 1 ), 2 ) ;
        xj = gNode( gElement( ie, 2 ), 1 ) ;
        yj = gNode( gElement( ie, 2 ), 2 ) ;
        xm = gNode( gElement( ie, 3 ), 1 ) ;
        ym = gNode( gElement( ie, 3 ), 2 ) ;
        ai = xj*ym - xm*yj ;
        aj = xm*yi - xi*ym ;
        am = xi*yj - xj*yi ;
        bi = yj - ym ;
        bj = ym - yi ;
        bm = yi - yj ;
        ci = -(xj-xm) ;
        cj = -(xm-xi) ;
        cm = -(xi-xj) ;
        area = (ai+aj+am)/2 ;
        B = [bi  0 bj  0 bm  0
              0 ci  0 cj  0 cm
             ci bi cj bj cm bm] ;
        B = B/2/area ;
    return

    function area = ElementArea( ie )
    %  计算单元面积
    %  输入参数:
    %     ie ----  单元号
    %  返回值:
    %     area  ----  单元面积
        global gNode gElement gMaterial

        xi = gNode( gElement( ie, 1 ), 1 ) ;
        yi = gNode( gElement( ie, 1 ), 2 ) ;
        xj = gNode( gElement( ie, 2 ), 1 ) ;
        yj = gNode( gElement( ie, 2 ), 2 ) ;
        xm = gNode( gElement( ie, 3 ), 1 ) ;
        ym = gNode( gElement( ie, 3 ), 2 ) ;
        ai = xj*ym - xm*yj ;
        aj = xm*yi - xi*ym ;
        am = xi*yj - xj*yi ;
        area = abs((ai+aj+am)/2) ;
    return

    function AssembleStiffnessMatrix( ie, k )
    %  把单元刚度矩阵集成到整体刚度矩阵
    %  输入参数:
    %      ie  --- 单元号
    %      k   --- 单元刚度矩阵
    %  返回值:
    %      无
        global gElement gK
        for i=1:1:3
            for j=1:1:3
                for p=1:1:2
                    for q=1:1:2
                        m = (i-1)*2+p ;
                        n = (j-1)*2+q ;
                        M = (gElement(ie,i)-1)*2+p ;
                        N = (gElement(ie,j)-1)*2+q ;
                        gK(M,N) = gK(M,N) + k(m,n) ;
                    end
                end
            end
        end
    return

    function egf = EquivalentGravityForce( ie )
    %  计算单元自重的等效节点力
    %  输入参数
    %      ie  ----- 单元号
    %  返回值
    %      egf ----- 自重的等效节点力
        global gElement gMaterial
       
        area = ElementArea( ie ) ;
        h    = gMaterial( gElement( ie, 4 ), 3 ) ;
        ro   = gMaterial( gElement( ie, 4 ), 4 ) ;
        w    = area * h * ro ;
        egf  = -w/3 * [0; 1; 0; 1; 0; 1] ;
    return

    function es = ElementStress( ie )
    %  计算单元的应力分量
    %  输入参数
    %      ie  ----- 单元号
    %  返回值
    %      es ----- 单元应力分量列阵(1×6): [sx, sy, txy, s1, s2, tmax]
        global gElement gDelta  gMaterial
       
        es = zeros( 1, 6 ) ;   % 单元的应力分量
        de = zeros( 6, 1 ) ;   % 单元节点位移列阵
        E  = gMaterial( gElement(ie, 4), 1 ) ;
        mu = gMaterial( gElement(ie, 4), 2 ) ;
        D = [ 1-mu    mu      0
               mu    1-mu     0
                0      0   (1-2*mu)/2] ;
        D = D*E/(1-2*mu)/(1+mu) ;
        B = MatrixB( ie ) ;
        for j=1:1:3
            de( 2*j-1 ) = gDelta( 2*gElement( ie, j )-1 ) ;
            de( 2*j   ) = gDelta( 2*gElement( ie, j )   ) ;
        end
        es(1:3) = D * B * de ;
        es(6) = 0.5*sqrt((es(1)-es(2))^2 + 4*es(3)^2 ) ;
        es(4) = 0.5*(es(1)+es(2)) + es(6) ;
        es(5) = 0.5*(es(1)+es(2)) - es(6) ;
    return

    function SaveResults( file_out )
    %  显示计算结果
    %  输入参数:
    %     file_out  --- 保存结果文件
    %  返回值:
    %     无

        global gNode gElement gMaterial gBC1 gDelta gNodeStress gElementStress
        save( file_out, 'gNode', 'gElement', 'gMaterial', 'gBC1', 'gDelta', 'gNodeStress', 'gElementStress' ) ;
    return

     

    展开全文
  • 在二维方向上的对圆柱体在空气域的流动问题,进行的网格划分程序
  • MATLAB编写三角形网格有限元程序,验证圣维南原理.
  • matlab naca0012网格划分

    2020-03-27 10:00:37
    matlab对naca0012翼型的网格划分文件,其中interfunction是翼型形状函数
  • Matlab 划分网格,三维曲面,网格划分好之后还需要提取网格上的每个点的位置坐标,之后需要用这些点的坐标进行运算。望各位大佬解惑,感激不尽

    Matlab 划分网格,三维曲面,网格划分好之后还需要提取网格上的每个点的位置坐标,之后需要用这些点的坐标进行运算。望各位大佬解惑,感激不尽

    展开全文
  • 输入ply点云进行三角网格matlab 输出三角网格化的效果
  • 将多边形分成由网格设置的更小的多边形 函数 PXY=DIVIDEXY(polygon,NX,NY) 输入多边形:由polygon.x(x坐标向量)和p​​olygon.y(y坐标向量)组成的结构NX:x 方向的分割数NY:y 方向的分割数 输出PXY:一个元胞...
  • MATLAB对三角网格进行线性细分

    千次阅读 2017-03-01 21:12:44
    这是除去三角网格的LOOP细分外, Jesus Mena的另一篇代码,是对三角网格进行线性细分,先贴代码 function [newVertices, newFaces] = linearSubdivision(vertices, faces) % Linear subdivision for triangle ...

    这是除去三角网格的LOOP细分外, Jesus Mena的另一篇代码,是对三角网格进行线性细分,先贴代码

    function [newVertices, newFaces] =  linearSubdivision(vertices, faces)
    % Linear subdivision for triangle meshes
    %
    %  Dimensions:
    %    vertices: 3xnVertices
    %    faces:    3xnFaces
    %  
    %  Author: Jesus Mena
    
    	global edgeVertex;
        global newIndexOfVertices;
    	newFaces = [];
    	newVertices = vertices;
    
    	nVertices = size(vertices,2);
    	nFaces    = size(faces,2);
    	edgeVertex= zeros(nVertices, nVertices);
    	newIndexOfVertices = nVertices;
    
        % ------------------------------------------------------------------------ %
    	% create a matrix of edge-vertices and a new triangulation (newFaces).
        % 
        % * edgeVertex(x,y): index of the new vertex between (x,y)
        %
        %  0riginal vertices: va, vb, vc.
        %  New vertices: vp, vq, vr.
        %
        %      vb                   vb             
        %     / \                  /  \ 
        %    /   \                vp--vq
        %   /     \              / \  / \
        % va ----- vc   ->     va-- vr --vc 
    	%
        
    	for i=1:nFaces
    		[vaIndex, vbIndex, vcIndex] = deal(faces(1,i), faces(2,i), faces(3,i));
    		
    		vpIndex = addEdgeVertex(vaIndex, vbIndex);
    		vqIndex = addEdgeVertex(vbIndex, vcIndex);
    		vrIndex = addEdgeVertex(vaIndex, vcIndex);
    		
    		fourFaces = [vaIndex,vpIndex,vrIndex; vpIndex,vbIndex,vqIndex; vrIndex,vqIndex,vcIndex; vrIndex,vpIndex,vqIndex]';
    		newFaces  = [newFaces, fourFaces]; 
        end;
        	
        % ------------------------------------------------------------------------ %
    	% positions of the new vertices
    	for v1=1:nVertices-1
    		for v2=v1:nVertices
    			vNIndex = edgeVertex(v1,v2);
                if (vNIndex~=0)
     				newVertices(:,vNIndex) = 1/2*(vertices(:,v1)+vertices(:,v2));
                end;
            end;
        end;
     	
    end
    
    % ---------------------------------------------------------------------------- %
    function vNIndex = addEdgeVertex(v1Index, v2Index)
    	global edgeVertex;
    	global newIndexOfVertices;
    
    	if (v1Index>v2Index) % setting: v1 <= v2
    		vTmp = v1Index;
    		v1Index = v2Index;
    		v2Index = vTmp;
    	end;
    	
    	if (edgeVertex(v1Index, v2Index)==0)  % new vertex
    		newIndexOfVertices = newIndexOfVertices+1;
    		edgeVertex(v1Index, v2Index) = newIndexOfVertices;
    	end;
    
    	vNIndex = edgeVertex(v1Index, v2Index);
    
        return;
    end
    命名为linearSubdivision.m文件。其代码的大致思想和之前LOOP细分是一样的,但是代码显然简单了不少,因为线性细分的规则是只需要在两点之间线性插值其中点,因此所用的代码简单非常多,相信不难理解,下面是用来绘图的plotMesh.m文件

    function plotMesh(vertices, faces)
        hold on;
        trimesh(faces', vertices(1,:), vertices(2,:), vertices(3,:));
        colormap hot;
        axis tight;
        axis square; 
        axis off;
        view(3);
    end
    最后是用来测试的test.m文件
     vertices = [10 10 10; -100 10 -10; -100 -10 10; 10 -10 -10]';
    faces = [1 2 3; 1 3 4; 1 4 2; 4 3 2]';
    
    figure(2);
    subplot(1,4,1);
    plotMesh(vertices, faces);
    for i=2:4
     subplot(1,4,i);
     [vertices, faces] = linearSubdivision(vertices, faces); 
     plotMesh(vertices, faces);
    end
    下面附上效果图

    线性细分四面体LOOP细分正四面体

    对比一下之前的LOOP细分,可以看到二者有明显的不同。

    展开全文
  • % mesh_sph 对给定的球壳进行网格划分,该球壳定义为% 半径 theta 方位角和 phi 极角。 欲了解更多信息%输入doc sph2cart。 Theta 是 0<theta<2*pi 和 pi/2<phi<pi/2 % rho 总是正双。 为了确定网格...
  • Delaunay算法的matlab实现,一种经典的三角网格划分方法。
  • 四面体(MATLAB)

    2017-06-09 21:16:01
    四面体网格划分的有限元计算源代码,在将变形体划分为四面体后,输入节点坐标,程序利用弹塑性本构,将载荷和位移关系转变为线性方程组,是一个很好的求解简单三维变形体的程序,还可输出单元变形量,应力等。
  • 三角形网格matlab有限元,matlab有限元网格划分,matlab源码.zip
  • 网格划分在很多领域都会使用,本次问题是插值需要采样,但区域是一个梯形。本文在此基础上进行了泛化,相关结论应该在CFD网格划分会有提及。本文应该只是重复造轮子。 问题描述 从矩形域的等距网格划分形变为任意...
  • Matlab二维插值加划分网格

    千次阅读 2019-08-22 16:49:53
    在均匀网格上插入散点数据 在均匀的查询点网格上插入随机分布的散点数据。 对函数介于 -2.5 和 2.5 之间的 200 个随机点采样。 xy = -2.5 + 5*gallery('uniformdata',[200 2],0); x = xy(:,1); y = xy(:,2); v = x.*...
  • r-square代码网格划分工具箱——meshpart 该工具箱包含用于多种图形和网格划分方法的 Matlab 代码,包括几何、光谱、几何光谱和坐标二分。 它还具有生成递归多路分区、顶点分隔符和嵌套解剖顺序的例程; 它有一些...
  • 该程序是基于邻域网格划分实现的聚类。版权原因,该程序内没有附带相应文章,大家可以根据代码内给定的标题和DOI号自行谷歌下载。 如果希望讨论聚类问题,欢迎根据代码内联系方式联系我,共同学习。如果只是初级的...
  • 可以对矩形、正方形、圆形板进行网格划分,获得的坐标和节点连接矩阵可用于进一步的有限元分析。 输入必要的输入,如长度、宽度、半径、角度、元素数量等,绘制有限元网格并获得节点连接性和节点坐标。 在板弯曲 ( ...
  • DistMesh, A Simple Mesh Generator in MATLAB, 有限元网格剖分不错工具,可做项目测试用,二维三维都可剖分,其主页介绍:http://persson.berkeley.edu/distmesh/
  • 1.语法:  [X,Y] = meshgrid(x,y)  [X,Y] = meshgrid(x)... 1.[X,Y] = meshgrid(x,y) 基于向量 x 和 y 中包含的坐标返回二维网格坐标。  X 是一个矩阵,每一行是 x 的一个副本;  Y 也是一个矩阵,每一列...
  • 只需下载文件,运行它,一个美丽的表面就在你面前。 您可以通过仅更改代码开头给出的X,Y,Z矩阵的元素来根据自己的需要更改表面,但切记不要更改这些元素的尺寸。 它最初必须是 4*4。 :-)

空空如也

空空如也

1 2 3 4 5 ... 20
收藏数 2,009
精华内容 803
关键字:

matlab网格划分

matlab 订阅