95 lines
2.6 KiB
Mathematica
95 lines
2.6 KiB
Mathematica
|
function [] = Tst()
|
||
|
|
% set up grid
|
||
|
|
gd=0.;
|
||
|
|
numElem=4;
|
||
|
|
eLen=0.25;
|
||
|
|
for i=1:numElem+1
|
||
|
|
ndCrd(i)=gd;
|
||
|
|
gd=gd+eLen;
|
||
|
|
end
|
||
|
|
for i=1:numElem
|
||
|
|
elemNod(i,1)=i;
|
||
|
|
elemNod(i,2)=i+1;
|
||
|
|
end
|
||
|
|
% initial level set
|
||
|
|
dpos=0.1;
|
||
|
|
for i=1:numElem+1
|
||
|
|
lSet(i)=sign(ndCrd(i)-dpos)*abs(dpos-ndCrd(i));
|
||
|
|
end
|
||
|
|
% force vector BC
|
||
|
|
fVect=zeros(12,1);
|
||
|
|
for i=1:numElem
|
||
|
|
if sign(lSet(elemNod(i,1)))~=sign(lSet(elemNod(i,2)))
|
||
|
|
fVect(elemNod(i,1))= 1.;
|
||
|
|
fVect(elemNod(i,2))= 1.;
|
||
|
|
end
|
||
|
|
end
|
||
|
|
% solve for velocity vector
|
||
|
|
Af=zeros(numElem+1);
|
||
|
|
% Assemble 'Stiffness' Matrices
|
||
|
|
for i=1:numElem
|
||
|
|
pos(1)=-1/sqrt(3);
|
||
|
|
pos(2)=1/sqrt(3);
|
||
|
|
AfL=zeros(2);
|
||
|
|
AfLGLS=zeros(2);
|
||
|
|
for j=1:2
|
||
|
|
shp(1)=(1-pos(j))/2.;
|
||
|
|
shp(2)=(1+pos(j))/2.;
|
||
|
|
dshp(1)=-0.5;
|
||
|
|
dshp(2)=0.5;
|
||
|
|
rset=shp(1)*lSet(elemNod(i,1))+shp(2)*lSet(elemNod(i,2));
|
||
|
|
dls=dshp(1)*lSet(elemNod(i,1))+dshp(2)*lSet(elemNod(i,2));
|
||
|
|
AfL=AfL+shp'*sign(rset)*(dls*dshp);
|
||
|
|
AfLGLS=AfLGLS+(dshp'*dls)*(1./abs(dls))*(dls*dshp);
|
||
|
|
end
|
||
|
|
for k=1:2;
|
||
|
|
for j=1:2;
|
||
|
|
Af(elemNod(i,j),elemNod(i,k))=Af(elemNod(i,j),elemNod(i,k))+AfL(j,k)+AfLGLS(j,k);
|
||
|
|
end
|
||
|
|
end
|
||
|
|
end
|
||
|
|
|
||
|
|
% Update level set
|
||
|
|
mMat=zeros(numElem+1);
|
||
|
|
mMatGLS=zeros(numElem+1);
|
||
|
|
f1=zeros(numElem+1,1);
|
||
|
|
f2=zeros(numElem+1,1);
|
||
|
|
f3=zeros(numElem+1,1);
|
||
|
|
for i=1:numElem
|
||
|
|
pos(1)=-1/sqrt(3);
|
||
|
|
pos(2)=1/sqrt(3);
|
||
|
|
mMatL=zeros(2);
|
||
|
|
mMatGLSL=zeros(2);
|
||
|
|
f1L=zeros(2,1);
|
||
|
|
f2L=zeros(2,1);
|
||
|
|
f3L=zeors(2,1);
|
||
|
|
for j=1:2
|
||
|
|
shp(1)=(1-pos(j))/2.;
|
||
|
|
shp(2)=(1+pos(j))/2.;
|
||
|
|
dshp(1)=-0.5;
|
||
|
|
dshp(2)=0.5;
|
||
|
|
Floc=shp(1)*fVect(elemNod(i,1))+shp(2)*fVect(elemNod(i,2));
|
||
|
|
rset=shp(1)*lSet(elemNod(i,1))+shp(2)*lSet(elemNod(i,2));
|
||
|
|
dls=dshp(1)*lSet(elemNod(i,1))+dshp(2)*lSet(elemNod(i,2));
|
||
|
|
mMatL=mMatL+shp'*shp;
|
||
|
|
mMatGLS=mMatGLS+dshp'*(dls/abs(dls))*Floc*1.*shp;
|
||
|
|
f1L=f1L+shp'*Floc*abs(dls);
|
||
|
|
f2L=f2L+dshp'*(dls/abs(dls))*Floc*1.*Floc*abs(dls);
|
||
|
|
vs=1.*((abs(0.1+Floc*abs(dls)))/(abs(Floc*abs(dls))+1.));
|
||
|
|
f3L=f3L+vs*dshp'*dls;
|
||
|
|
end
|
||
|
|
for k=1:2;
|
||
|
|
for j=1:2;
|
||
|
|
mMat(elemNod(i,j),elemNod(i,k))=mMat(elemNod(i,j),elemNod(i,k))+mMatL(j,k);
|
||
|
|
mMatGLS(elemNod(i,j),elemNod(i,k))=mMatGLS(elemNod(i,j),elemNod(i,k))+mMatGLSL(j,k);
|
||
|
|
end
|
||
|
|
f1(elemNod(i,k))=f1(elemNod(i,k))+f1L(k);
|
||
|
|
f2(elemNod(i,k))=f2(elemNod(i,k))+f2L(k);
|
||
|
|
f3(elemNod(i,k))=f3(elemNod(i,k))+f3L(k);
|
||
|
|
end
|
||
|
|
end
|
||
|
|
dt=0.1;
|
||
|
|
lSet=lSet-(((mMat+mMatGLS)^-1)/dt)*(f1+f2+f3);
|
||
|
|
lSet
|
||
|
|
|
||
|
|
|