Mesh_remi.zip

function [coord edof dofs meshdb]=delaunay2d(vertices,segments,mp,holes) % [coord edof dofs]=delaunay2d(vertices,segments,mp,holes) %------------------------------------------------------------- % PURPOSE % calculates the delaunay mesh by using randomized incremental insertion % triangulation on input followed by adding steiner point until the area % criteria is fullfilled and/or bisection split of triangles of low % quality. % % INPUT: segments : boundary topology based on segments, nseg x 2 % vertices : the global coordinates of vertices, nv x 2 % holes : optional: struct of segments of holes, {nseg x 2} % mp : unstructured mesh [ maxArea dpn nen fixed qmin], % % nel : number of index elements % ned : number of element dofs % np : number of nodes in mesh % % OUTPUT: edof : topology matrix , nel x ned+1 % coord : global coordinate matrix, np x 2 % dofs : global nodal dof matrix, np x dpn %------------------------------------------------------------- % DEPENDENCY: jacobi,changedofs %------------------------------------------------------------- % LAST MODIFIED: J.Lorentzon, Jonas Lindemann, Kent Persson 2010-04-18 % Copyright (c) Division of Structural Mechanics and % Department of Solid Mechanics. % Lund Institute of Technology %------------------------------------------------------------- global T;global tri;global eps; % if nargin<4;holes=[];end; eps=1e-12; % NOTE: tri is a item holder, T is an cellarray tri=[];T={}; % delaunay randomized incremental algorithm triangulation Triangulate(vertices); % remove outside triangle using ray trace technique trim(segments,holes); % set convergence criteria maxArea=mp(1)*(mp(3)-2); % fixed boundary? if numel(mp)>=4 addbp=mp(4); else addbp=1; end % limit split only to given quality factor q of elements qlimit=0.5;Jm=2*maxArea; if numel(mp)==5 qmin=mp(5); else qmin=0.6; end % iterative procedure steiner=0;bisect=0;ntotal=0;enforce=0;enter=0; while maxArea<max(abs(Jm)) Jm=0;counter=enter+bisect; for i=1:numel(T) if isempty(T{i}.child) && T{i}.level>0 ntotal=ntotal+1; [q J]=qfactor(T{i}.edof); Jm=max(Jm,abs(J)); parent=i;Pi=T{i}.edof(1);Pj=T{i}.edof(2);Pl=T{i}.edof(3); if q>qlimit && maxArea<abs(J) && q rel='nofollow' onclick='return false;'>0 || (maxArea<abs(J) && q rel='nofollow' onclick='return false;'>0 && enforce) enforce=0;enter=enter+1; % for regular triangle, Steiner point is prefered! nel=numel(T); tri.vertices(end+1,:)=mean(tri.vertices([Pi Pj Pl],:)); Pr=size(tri.vertices,1); % T{parent}.child=[ nel+1 nel+2 nel+3]; edge=OppositeEdge([parent parent parent],[Pl Pi Pj]); % create the new three triangles T{nel+1}=CreateTriangle([Pr,Pi,Pj],nel,parent,nel+1,[nel+3 T{parent}.e(edge(1)) nel+2],[Pl T{parent}.o(edge(1)) Pl]); T{nel+2}=CreateTriangle([Pr,Pj,Pl],nel,parent,nel+2,[nel+1 T{parent}.e(edge(2)) nel+3],[Pi T{parent}.o(edge(2)) Pi]); T{nel+3}=CreateTriangle([Pr,Pl,Pi],nel,parent,nel+3,[nel+2 T{parent}.e(edge(3)) nel+1],[Pj T{parent}.o(edge(3)) Pj]); % the insertion affect the enviroments adjacent list as well. for j=1:3 % ignore null elements (e=0) if T{parent}.e(edge(j)) [dummy, i1]=find(T{T{parent}.e(edge(j))}.e==parent); T{T{parent}.e(edge(j))}.e(i1)=nel+j;T{T{parent}.e(edge(j))}.o(i1)=Pr; end; end; % enforce delaunay condition LegalizeEdge(T{nel+1}); LegalizeEdge(T{nel+2}); LegalizeEdge(T{nel+3}); steiner=steiner+1; elseif q<qlimit && maxArea<abs(J) && q rel='nofollow' onclick='return false;'>0 % bisection, purpose, improve the q factor followed by % flipping edges I=T{i}.edof; indp=[2 3 1];indm=[3 1 2]; edges=[ I(1) I(2); I(2) I(3);I(3) I(1)]; % identify longest edge and add midpoint [dummy, i4]=max(sum((tri.vertices(edges(:,1),:)-tri.vertices(edges(:,2),:)).^2,2)); if i4==1; % along edge 1 P1=Pi;P2=Pj;P3=Pl;m=1; elseif i4==3; % along edge 3 P1=Pl;P2=Pi;P3=Pj;m=3; else % along edge 2 P1=Pj;P2=Pl;P3=Pi;m=2; end; % identify enviroment ek=T{parent}.e(m);Pk=T{parent}.o(m); edge=OppositeEdge(T{parent}.e(m),T{parent}.o(m));e1=0;e2=0;o1=0;o2=0; if Pk<=0; enforce=1;end; if ek;f=T{ek}.level>0;else f=0;enforce=1;end; if (~ek && addbp) || ek && ((Pk>0 && f)||addbp) enter=enter+1; tri.vertices(end+1,:)=mean(tri.vertices([P1 P2],:)); Pr=size(tri.vertices,1); if ek ; e1=T{ek}.e(indp(edge));e2=T{ek}.e(indm(edge)); end; if ek ; o1=T{ek}.o(indp(edge));o2=T{ek}.o(indm(edge)); end; % child update nel=numel(T);fixed=nel; if ek ; if T{ek}.level<=0; fixed=0;end;end; T{parent}.child=[ nel+1 nel+2]; if ek ; T{ek}.child=[ nel+3 nel+4];end; % build the four new triangles T{nel+1}=CreateTriangle([ Pr P2 P3],nel,parent,nel+1,[(nel+4)*(~~ek) T{parent}.e(indp(m)) nel+2],[Pk*(~~ek) T{parent}.o(indp(m)) P1]); T{nel+2}=CreateTriangle([ Pr P3 P1],nel,parent,nel+2,[nel+1 T{parent}.e(indm(m)) (nel+3)*(~~ek)],[P2 T{parent}.o(indm(m)) Pk*(~~ek) ]); if ek; T{nel+3}=CreateTriangle([ Pr P1 Pk],fixed,ek,nel+3,[nel+2 e1 nel+4],[P3 o1 P2]);end; if ek; T{nel+4}=CreateTriangle([ Pr Pk P2],fixed,ek,nel+4,[nel+3 e2 nel+1],[P1 o2 P3]);end; % the insertion affect the enviroments adjacent list as well. % ek element enviroment edge=OppositeEdge([e1 e2],[o1 o2]); if e1 ; T{e1}.o(edge(1))=Pr;T{e1}.e(edge(1))=nel+3;end; if e2 ; T{e2}.e(edge(2))=nel+4;T{e2}.o(edge(2))=Pr;end; % parent enviroment e3=T{parent}.e(indp(m));e4=T{parent}.e(indm(m)); edge=OppositeEdge([e3 e4],[T{parent}.o(indp(m)) T{parent}.o(indm(m))]); if e3 ; T{e3}.e(edge(1))=nel+1;T{e3}.o(edge(1))=Pr;end; if e4 ; T{e4}.e(edge(2))=nel+2;T{e4}.o(edge(2))=Pr;end; % enforce delaunay condition LegalizeEdge(T{nel+1}); LegalizeEdge(T{nel+2}); if ek LegalizeEdge(T{nel+3}); LegalizeEdge(T{nel+4}); end; bisect=bisect+1; end; % end end end qlimit=min(max((steiner+1)/(bisect+1)*qlimit*(enter/ntotal),qmin),0.9); if (enter+bisect)==counter Jm=0; end end % convert into CALFEM format [coord edof dofs]=extractT; % set mesh property if mp(3)==4 [coord edof dofs]= convertMesh(edof,coord,dofs); end % set dofspernode [edof dofs ]=changedofs(edof,dofs, mp(2)); % clear globals! meshdb=T; % T=[];tri=[]; % finished! function [q J]=qfactor(etop) % q factor is the quotient product of the angle with optimal angle % the area J of triangle is a byproduct suited to current algorithm. global tri; % if min(etop)<1 q=-1; J=1