restart;with(plots):
setoptions3d(title=`Woven Unit Cell`,style=PATCH, axes=BOXED);
H:=0.00635: L:=0.015: W:=0.015:
NH:=12: NL:=16: NW:=16:
a:=0.0028: b:=0.0009: c:=0.001:
f1p:=[L/4+a*sin(t),y,H/2+c*sin(2*Pi*y/W)+b*cos(t)]:
f2p:=[3*L/4+a*sin(t),y,H/2+c*sin(2*Pi*y/W+Pi)+b*cos(t)]:
f3p:=[x,3*W/4+a*sin(t),H/2+c*sin(2*Pi*x/L)+b*cos(t)]:
f4p:=[x,W/4+a*sin(t),H/2+c*sin(2*Pi*x/L+Pi)+b*cos(t)]:
fiber1:=plot3d([f1p[1],f1p[2],f1p[3]],y=0..W,t=0..2*Pi,color=green):
fiber2:=plot3d([f2p[1],f2p[2],f2p[3]],y=0..W,t=0..2*Pi,color=green):
fiber3:=plot3d([f3p[1],f3p[2],f3p[3]],x=0..L,t=0..2*Pi,color=yellow):
fiber4:=plot3d([f4p[1],f4p[2],f4p[3]],x=0..L,t=0..2*Pi,color=yellow):
display(fiber1,fiber2,fiber3,fiber4);
f1ineq:=((x-L/4)/a)^2+((z-(H/2+c*sin(2*Pi*y/W)))/b)^2-1:
f2ineq:=((x-3*L/4)/a)^2+((z-(H/2+c*sin(2*Pi*y/W+Pi)))/b)^2-1:
f3ineq:=((y-3*W/4)/a)^2+((z-(H/2+c*sin(2*Pi*x/L)))/b)^2-1:
f4ineq:=((y-W/4)/a)^2+((z-(H/2+c*sin(2*Pi*x/L+Pi)))/b)^2-1:
elements:=array(1..NL*NW*NH):
for E from 1 to NL*NW*NH do elements[E]:=0 od:
ELEMENT:=array(1..NL,1..NW,1..NH):
for i from 1 to NL do
for j from 1 to NW do
for k from 1 to NH do
ELEMENT[i,j,k]:=0
od:od:od:
nodes:=array(1..NL*NW*NH):
X:=array(1..(NL+1),1..(NW+1),1..(NH+1)):
for i from 1 to NL+1 do
for j from 1 to NW+1 do
for k from 1 to NH+1 do
X[i,j,k]:=(i-1)*L/NL:
od:od:od:
Y:=array(1..(NL+1),1..(NW+1),1..(NH+1)):
for i from 1 to NL+1 do
for j from 1 to NW+1 do
for k from 1 to NH+1 do
Y[i,j,k]:=(j-1)*W/NW:
od:od:od:
Z:=array(1..(NL+1),1..(NW+1),1..(NH+1)):
for i from 1 to NL+1 do
for j from 1 to NW+1 do
for k from 1 to NH+1 do
Z[i,j,k]:=(k-1)*H/NH:
od:od:od:
for i from 1 to NL do
for j from 1 to NW do
for k from 1 to NH do
EL_NUM:=i+(j-1)*NL+(k-1)*NW*NL;
nodes[EL_NUM]:=[i+(j-1)*(NL+1)+(k-1)*(NW+1)*(NL+1),
(i+1)+(j-1)*(NL+1)+(k-1)*(NW+1)*(NL+1),
(i+1)+j*(NL+1)+(k-1)*(NW+1)*(NL+1),
i+j*(NL+1)+(k-1)*(NW+1)*(NL+1),
i+(j-1)*(NL+1)+k*(NW+1)*(NL+1),
(i+1)+(j-1)*(NL+1)+k*(NW+1)*(NL+1),
(i+1)+j*(NL+1)+k*(NW+1)*(NL+1),
i+j*(NL+1)+k*(NW+1)*(NL+1)];
XCEN:=(i-1/2)*L/NL;
YCEN:=(j-1/2)*W/NW;
ZCEN:=(k-1/2)*H/NH;
TOL:=((L*H*W)^(1/3))/(NL*NW*NH);
if(evalf(subs(x=XCEN,y=YCEN,z=ZCEN,f1ineq))<TOL) then
elements[EL_NUM]:=1 fi;
if(evalf(subs(x=XCEN,y=YCEN,z=ZCEN,f2ineq))<TOL) then
elements[EL_NUM]:=2 fi;
if(evalf(subs(x=XCEN,y=YCEN,z=ZCEN,f3ineq))<TOL) then
elements[EL_NUM]:=3 fi;
if(evalf(subs(x=XCEN,y=YCEN,z=ZCEN,f4ineq))<TOL) then
elements[EL_NUM]:=4 fi;
od: od: od:
SUM:=0; for E from 1 to NL*NW*NH do
if (elements[E]=0) then SUM:=evalf (SUM)+1;
fi;
od;
Vf:=1-SUM/NL/NW/NH;
Nnum:=(NL+1)*(NW+1)*(NH+1);
Enum:=NL*NW*NH;
for i from 1 to NL do
for j from 1 to NW do
for k from 1 to NH do
NINTS:=0;
I1:=i; I2:=i+1; I3:=i+1; I4:=i; J1:=j; J2:=j; J3:=j+1; J4:=j+1; K1:=k;
K2:=k; K3:=k; K4:=k;
I5:=i; I6:=i+1; I7:=i+1; I8:=i; J5:=j; J6:=j; J7:=j+1; J8:=j+1;
K5:=k+1; K6:=k+1; K7:=k+1; K8:=k+1;
if(evalf(subs(x=X[I1,J1,K1],y=Y[I1,J1,K1],z=Z[I1,J1,K1],f1ineq))<TOL)
then
NINTS:=NINTS+1; fi:
if(evalf(subs(x=X[I2,J2,K2],y=Y[I2,J2,K2],z=Z[I2,J2,K2],f1ineq))<TOL)
then
NINTS:=NINTS+1; fi:
if(evalf(subs(x=X[I3,J3,K3],y=Y[I3,J3,K3],z=Z[I3,J3,K3],f1ineq))<TOL)
then
NINTS:=NINTS+1; fi:
if(evalf(subs(x=X[I4,J4,K4],y=Y[I4,J4,K4],z=Z[I4,J4,K4],f1ineq))<TOL)
then
NINTS:=NINTS+1; fi:
if(evalf(subs(x=X[I5,J5,K5],y=Y[I5,J5,K5],z=Z[I5,J5,K5],f1ineq))<TOL)
then
NINTS:=NINTS+1; fi:
if(evalf(subs(x=X[I6,J6,K6],y=Y[I6,J6,K6],z=Z[I6,J6,K6],f1ineq))<TOL)
then
NINTS:=NINTS+1; fi:
if(evalf(subs(x=X[I7,J7,K7],y=Y[I7,J7,K7],z=Z[I7,J7,K7],f1ineq))<TOL)
then
NINTS:=NINTS+1; fi:
if(evalf(subs(x=X[I8,J8,K8],y=Y[I8,J8,K8],z=Z[I8,J8,K8],f1ineq))<TOL)
then
NINTS:=NINTS+1; fi:
ELEMENT1[i,j,k]:=NINTS;
od;od;od;
for i from 1 to NL do
for j from 1 to NW do
for k from 1 to NH do
NINTS:=0;
I1:=i; I2:=i+1; I3:=i+1; I4:=i; J1:=j; J2:=j; J3:=j+1; J4:=j+1;
K1:=k; K2:=k; K3:=k; K4:=k;
I5:=i; I6:=i+1; I7:=i+1; I8:=i; J5:=j; J6:=j; J7:=j+1; J8:=j+1;
K5:=k+1; K6:=k+1; K7:=k+1; K8:=k+1;
if(evalf(subs(x=X[I1,J1,K1],y=Y[I1,J1,K1],z=Z[I1,J1,K1],f2ineq))<TOL)
then
NINTS:=NINTS+1; fi:
if(evalf(subs(x=X[I2,J2,K2],y=Y[I2,J2,K2],z=Z[I2,J2,K2],f2ineq))<TOL)
then
NINTS:=NINTS+1; fi:
if(evalf(subs(x=X[I3,J3,K3],y=Y[I3,J3,K3],z=Z[I3,J3,K3],f2ineq))<TOL)
then
NINTS:=NINTS+1; fi:
if(evalf(subs(x=X[I4,J4,K4],y=Y[I4,J4,K4],z=Z[I4,J4,K4],f2ineq))<TOL)
then
NINTS:=NINTS+1; fi:
if(evalf(subs(x=X[I5,J5,K5],y=Y[I5,J5,K5],z=Z[I5,J5,K5],f2ineq))<TOL)
then
NINTS:=NINTS+1; fi:
if(evalf(subs(x=X[I6,J6,K6],y=Y[I6,J6,K6],z=Z[I6,J6,K6],f2ineq))<TOL)
then
NINTS:=NINTS+1; fi:
if(evalf(subs(x=X[I7,J7,K7],y=Y[I7,J7,K7],z=Z[I7,J7,K7],f2ineq))<TOL)
then
NINTS:=NINTS+1; fi:
if(evalf(subs(x=X[I8,J8,K8],y=Y[I8,J8,K8],z=Z[I8,J8,K8],f2ineq))<TOL)
then
NINTS:=NINTS+1; fi:
ELEMENT2[i,j,k]:=NINTS;
od;od;od;
for i from 1 to NL do
for j from 1 to NW do
for k from 1 to NH do
NINTS:=0;
I1:=i; I2:=i+1; I3:=i+1; I4:=i; J1:=j; J2:=j; J3:=j+1; J4:=j+1;
K1:=k; K2:=k; K3:=k; K4:=k;
I5:=i; I6:=i+1; I7:=i+1; I8:=i; J5:=j; J6:=j; J7:=j+1; J8:=j+1;
K5:=k+1; K6:=k+1; K7:=k+1; K8:=k+1;
if(evalf(subs(x=X[I1,J1,K1],y=Y[I1,J1,K1],z=Z[I1,J1,K1],f3ineq))<TOL)
then
NINTS:=NINTS+1; fi:
if(evalf(subs(x=X[I2,J2,K2],y=Y[I2,J2,K2],z=Z[I2,J2,K2],f3ineq))<TOL)
then
NINTS:=NINTS+1; fi:
if(evalf(subs(x=X[I3,J3,K3],y=Y[I3,J3,K3],z=Z[I3,J3,K3],f3ineq))<TOL)
then
NINTS:=NINTS+1; fi:
if(evalf(subs(x=X[I4,J4,K4],y=Y[I4,J4,K4],z=Z[I4,J4,K4],f3ineq))<TOL)
then
NINTS:=NINTS+1; fi:
if(evalf(subs(x=X[I5,J5,K5],y=Y[I5,J5,K5],z=Z[I5,J5,K5],f3ineq))<TOL)
then
NINTS:=NINTS+1; fi:
if(evalf(subs(x=X[I6,J6,K6],y=Y[I6,J6,K6],z=Z[I6,J6,K6],f3ineq))<TOL)
then
NINTS:=NINTS+1; fi:
if(evalf(subs(x=X[I7,J7,K7],y=Y[I7,J7,K7],z=Z[I7,J7,K7],f3ineq))<TOL)
then
NINTS:=NINTS+1; fi:
if(evalf(subs(x=X[I8,J8,K8],y=Y[I8,J8,K8],z=Z[I8,J8,K8],f3ineq))<TOL)
then
NINTS:=NINTS+1; fi:
ELEMENT3[i,j,k]:=NINTS;
od;od;od;
for i from 1 to NL do
for j from 1 to NW do
for k from 1 to NH do
NINTS:=0;
I1:=i; I2:=i+1; I3:=i+1; I4:=i; J1:=j; J2:=j; J3:=j+1; J4:=j+1;
K1:=k; K2:=k; K3:=k; K4:=k;
I5:=i; I6:=i+1; I7:=i+1; I8:=i; J5:=j; J6:=j; J7:=j+1; J8:=j+1;
K5:=k+1; K6:=k+1; K7:=k+1; K8:=k+1;
if(evalf(subs(x=X[I1,J1,K1],y=Y[I1,J1,K1],z=Z[I1,J1,K1],f4ineq))<TOL)
then
NINTS:=NINTS+1; fi:
if(evalf(subs(x=X[I2,J2,K2],y=Y[I2,J2,K2],z=Z[I2,J2,K2],f4ineq))<TOL)
then
NINTS:=NINTS+1; fi:
if(evalf(subs(x=X[I3,J3,K3],y=Y[I3,J3,K3],z=Z[I3,J3,K3],f4ineq))<TOL)
then
NINTS:=NINTS+1; fi:
if(evalf(subs(x=X[I4,J4,K4],y=Y[I4,J4,K4],z=Z[I4,J4,K4],f4ineq))<TOL)
then
NINTS:=NINTS+1; fi:
if(evalf(subs(x=X[I5,J5,K5],y=Y[I5,J5,K5],z=Z[I5,J5,K5],f4ineq))<TOL)
then
NINTS:=NINTS+1; fi:
if(evalf(subs(x=X[I6,J6,K6],y=Y[I6,J6,K6],z=Z[I6,J6,K6],f4ineq))<TOL)
then
NINTS:=NINTS+1; fi:
if(evalf(subs(x=X[I7,J7,K7],y=Y[I7,J7,K7],z=Z[I7,J7,K7],f4ineq))<TOL)
then
NINTS:=NINTS+1; fi:
if(evalf(subs(x=X[I8,J8,K8],y=Y[I8,J8,K8],z=Z[I8,J8,K8],f4ineq))<TOL)
then
NINTS:=NINTS+1; fi:
ELEMENT4[i,j,k]:=NINTS;
od;od;od;
minimum1:=proc(x,yp,z)
MIND:=10^100;
MMT=-100;
MMY:=-100;
M:=20;
for mm from 0 to M do
for nn from 0 to M do
tm:=(mm/M)*2*evalf(Pi);
ym:=(nn/M)*W;
dist:=evalf((x-subs(t=tm,f1p[1]))^2+
(yp-subs(y=ym,f1p[2]))^2+
(z-subs(y=ym,t=tm,f1p[3]))^2);
if(dist<MIND) then
MIND:=dist;
MMT:=tm;
MMY:=ym;
fi;
od;od;
RETURN (evalf(subs(t=MMT,f1p[1])),
evalf(subs(y=MMY,f1p[2])),
evalf(subs(y=MMY,t=MMT,f1p[3])));
end:
g1:=minimum1;
minimum2:=proc(x,yp,z)
MIND:=10^100;
MMT=-100;
MMY:=-100;
M:=20;
for mm from 0 to M do
for nn from 0 to M do
tm:=(mm/M)*2*evalf(Pi);
ym:=(nn/M)*L;
dist:=evalf((x-subs(t=tm,f2p[1]))^2+
(yp-subs(y=ym,f2p[2]))^2+
(z-subs(y=ym,t=tm,f2p[3]))^2);
if(dist<MIND) then
MIND:=dist;
MMT:=tm;
MMY:=ym;
fi;
od;od;
RETURN (evalf(subs(t=MMT,f2p[1])),
evalf(subs(y=MMY,f2p[2])),
evalf(subs(y=MMY,t=MMT,f2p[3])));
end:
g2:=minimum2;
minimum3:=proc(xp,y,z)
MIND:=10^100;
MMT=-100;
MMX:=-100;
M:=20;
for mm from 0 to M do
for nn from 0 to M do
tm:=(mm/M)*2*evalf(Pi);
xm:=(nn/M)*L;
dist:=evalf((xp-subs(x=xm,f3p[1]))^2+
(y-subs(t=tm,f3p[2]))^2+
(z-subs(x=xm,t=tm,f3p[3]))^2);
if(dist<MIND) then
MIND:=dist;
MMT:=tm;
MMX:=xm;
fi;
od;od;
RETURN (evalf(subs(x=MMX,f3p[1])),
evalf(subs(t=MMT,f3p[2])),
evalf(subs(x=MMX,t=MMT,f3p[3])));
end:
g3:=minimum3;
minimum4:=proc(xp,y,z)
MIND:=10^100;
MMT=-100;
MMX:=-100;
M:=20;
for mm from 0 to M do
for nn from 0 to M do
tm:=(mm/M)*2*evalf(Pi);
xm:=(nn/M)*L;
dist:=evalf((xp-subs(x=xm,f4p[1]))^2+
(y-subs(t=tm,f4p[2]))^2+
(z-subs(x=xm,t=tm,f4p[3]))^2);
if(dist<MIND) then
MIND:=dist;
MMT:=tm;
MMX:=xm;
fi;
od;od;
RETURN (evalf(subs(x=MMX,f4p[1])),
evalf(subs(t=MMT,f4p[2])),
evalf(subs(x=MMX,t=MMT,f4p[3])));
end:
g4:=minimum4;