Matlab网格划分程序Distmesh讲解.docx
Matlab网格划分程序Distmesh讲解（一）DiMmeSh是个maUab语言写的W格划分软件.源文件可以从上面的网址获褥.这电按行讲解各个邦例。PO1.dCm概宛例是一个雎位圈（半径为1）的网格划分,划分后的网格为：以下逐行讲解该比例：IunclionPo1.demo(itcration_max.h)%Parameters:%Input,integer11iRATI0N_MAX.themaximumnumberofiterationsthatD1STMESH%shouldtake.(Theprogrammightlakefeweriterationsifitdetcc<sconvergence.)%Input,realh.theneshspacingPanImeler.%这里行两个输入参数.一个是门ERATION_MAX.迭代的最大次数.%另一个是h.网格划分的大小，(Xh<l%院认参数值为：FrERATlON=200h=0.lp,l=distmesh_2d(fd,ih.h,box.iieralion_max.txed);函数须要至少六个参数.d=M(p).p=xWfa给定任一点到边界的距感函数,本例中定义为:<1=sqrt(xA2+yA2).|;fl),sealededgelengthfunctionh(x.y).也就是卜4J格大小的函数.H)也就是h,网格的大小realB()X(2.2).theboundingboxxmin.ymin;XmaX.ymax.最大外围矩形范附本例中为0.0J.IJITERATIONeMAX.themaximumnumberofiterations.realPFlX(NFIX,2Ihefixednodepositions.网格中须要固定的点坐标,也就是肯定须耍出现在网格中的点.输出参数：malP(N,2)Ihenodepositions.网格点的x.y坐标integerT(NT,3).thetriangleindices.输出网格任个二角形的.个原点.第步：x.y|=ncshgrid(box(l,l):h():box(2J),.box(12):hO-qrt(3V2:box(2.2)M褥hO,网格的大小,先把能涵盖欲划分区域的垃大矩形划分为结构网格.然后把偶数行的点整体向右平移半格.08其次步:依据W的函数定义，移除边界外的点。p=p(fcx,alj<fd.p.vrargi11:)<=geps,:);varaginfd.Ai的Pft加参数.这里为空.gcps=0(X)lh:也就是保留了到边界的距离以外O<K)l-h以内的点.密度函数fh,每个点上产生一个O随机数，推断是否小于rtVmax（八）)大于的话,改点被删除。p=(pllx:p(rand(sizc(p.I),1)<rt)./max(rt>),:);InHX.dummy=SiZC(Pfix):当指定f某些点要保留的时候,把保留的点加入,删除型复的点。%Especiallywhentheuserhasilud<xifixedpoints,wemayhaveafew%duplicates.Gclridofanyyouspot.%p=unique(p.,rows,);N=sizc(p.IX这个时候产生的网格如下：代PoId=inf;%第次迭代前设置旧的点的坐标为无力while(iteration<iteration_max)iteration=iteration1;%先推断上次移动后的点和旧的点之间的移动距离，假如小于某个阀位，悴止迭代if(UOl<max(Sqrt(sum(p-pold).2,2)h)pold=p:%假如还可以移动,保存当我的节点I=delaunayn(P):%利用delauny算法,生成三角形网格Iriangulatioiucount=(riangulaiion.)unt+Izpmid=(+p(tC2).:)+p<t(:.3).:)/3;%计算三角形的心t=(fcval_r(f,pmid,varargin(:)<=-gcps,:);%移除曳心在边界外部的珀形%4,Describeeachbarbyauniquepairofnodes.%生成网格的边的集合.也就是相邻点之间连接的线段bars=(:J23|)Jibars=unique(sort(bars.2).'rows*);end%6,Moven)eshpointsbasedonbarlengths1.andforcesF%Makealistofthebarvectorsandlengths.%Sel1.Otothedesiredlengths.FtotheSCHiIrbarforces.%andFVECtothex,yCOmPonCntSofthebarforces.%Allhcfixedpositions,resettheforce100.%barvcc-p(ban-p(bars(:.2),:);%生成bar的矢Kt1.=Sqft(sum(barvec?2.2);%itWbar的长度%依据每个bar的中点坐标il算须要的三角形边的边长(这个在fh函数里限制,hbrs=fcva1.tfh,(p(t>ars(:J),:)*p(bais(:,2),:)/2,varargin:);%计算须要的bar的长度，已经乘上了两个scale参数Fscale.sqrl(surn(1.2)/sum(hbm.A2):%详细可参考他们的PaPCr1.O=hbarstFscale*8qr(SUm(1.F2sum(hbais.A2);%计算斑个Mr上力F=max(1.O-1.,0);%bar上力的筮此x.y方向Fvec=F./1.11.11.*barvec;%计算Ftrn.每个节点上力的残量FWt=full(sparsc(hars(:JI.112,2),oncs(sizc(F)1,2,1.2J4Fvcc,-Fvcc,N,2);%对于固定点，力的残量为零Roi(ksizc(pfixJ),:)=0:%依据每个节点上的受力.移动该点p=+deha*Fo1:%7.BnngoutsidepointsbackIotheboundary%Use(henumericalgradientofFDoxojecPOimSback(oIheboundary.d=fcval_r(M.p.VarargiII:);%il算点到边界距国ix=d>0;%it算移动梯度,相对边界dgradx=(feval_r(fd4p<ix.l)+deps.p(ix.2).vanirgin|:|>-d(ix)/dcpsdgra<ly-(feval_r(fd,|p<ixj),p(ix,2)*dcps|,varargin|:)-d(ix)/dcps;%将这些扬动到边界外的投射回边界上p(ix,:)=p(ix.:)-1d(ix).dgradx.d(ix).4dgrady:%Ineeded(hefollowingrndifivationtof>rccthefixedpoints(oslay.%Othcnvisc.theycandriftoutsidetheregionandbelost.%JVB.21AugUM2008.%p(linfix.1:2)=pfix:N=size(p,I);%8.Terminationcriterion:.Allinteriorn<lcsmovelessthandp<ol(sealed)%if(max(sqrt(sum(de)tat*Ftol(d<-geps.:).A2.2)/h)<dplol)break:endend
