printf("\n bpatch.txt routines for creating and analyzing Bezier patches");
printf("\nwith(plots): with(plottools): with(linalg):");
with(plots): with(plottools): with(linalg):

printf("Commands in patch.txt");
printf("\n  BPatch(int m, int n, matrix pts) - computes a Bezier patch from a matrix of control points");
printf("\n  PatchPlot(vector patch) - plots a Bezier patch with an isoparametric net");
printf("\n  ContNet(int m, int n, matrix pts) - plots the control net of a bicubic Bezier patch");
printf("\n           from a matrix of control points");
printf("\n  InterPatch(int m, int n, vector U, vector V,  matrix datapts) - computes the control points for a Bezier patch");
printf("\n           that interpolates the matrix of pts in mxn grid");
printf("\n  LstSqrFit(int m, int n, vector U, vector V, matrix datapts, int M, int N) - computes the contol points for a Bezier patch");
printf("\n           that best fits the data points in pts");
printf("\n  SplinePatch(pts) - creates a C2 Bezier patch from a matrix of spline control points of form [i,j,x(i,j),y(i,j),z(i,j)]");
printf("\n  SplineNet(pts) - displays the control net for a C2 composite Bezier patch frpm spline control points");
printf("\n SplineNetA(pts) - displays the spline control net for a C2 composite Bezier patch");

BPatch:=proc (m,n,pts)
  local i,j,k,XCrd,YCrd,ZCrd,X,Y,Z,Bi,Bj,ans;

  XCrd:= array(0..m,0..n):
  YCrd:= array(0..m,0..n):
  ZCrd:= array(0..m,0..n):
  for i from 0 to m do
    for j from 0 to n do
      k:=1+i+(m+1)*j:
      XCrd[i,j]:=pts[k,1]:
      YCrd[i,j]:=pts[k,2]:
      ZCrd[i,j]:=pts[k,3]:
    od:
  od:
  unassign('i','j','k'):

  X:=simplify(sum(binomial(m,i)*u^i*(1-u)^(m-i)*sum(binomial(n,j)*v^j*(1-v)^(n-j)*XCrd[i,j],j=0..n),i=0..m)):
  Y:=simplify(sum(binomial(m,i)*u^i*(1-u)^(m-i)*sum(binomial(n,j)*v^j*(1-v)^(n-j)*YCrd[i,j],j=0..n),i=0..m)):
  Z:=simplify(sum(binomial(m,i)*u^i*(1-u)^(m-i)*sum(binomial(n,j)*v^j*(1-v)^(n-j)*ZCrd[i,j],j=0..n),i=0..m)):

  ans:=vector([X,Y,Z]):
  RETURN (evalm(ans))
end:
 
PatchPlot:=proc(patch)
   plot3d(patch,u=0..1,v=0..1,scaling=constrained,axes=boxed,tickmarks=[0,0,0]):
end:

ContNet:=proc(m,n,pts)
  local i,j,k,b,p,q,A,B,C:

  b:=array(0..m,0..n):
  for i from 0 to m do
    for j from 0 to n do
      k:=1+i+(m+1)*j:
      b[i,j]:=[pts[k,1],pts[k,2],pts[k,3]]:
    od:
  od:
  unassign('i','j'):
  for i from 0 to m do
    p[i]:=curve([seq(b[i,j],j=0..n)],color=black):
  od:
  unassign('i'):
  for j from 0 to n do
    q[j]:=curve([seq(b[i,j],i=0..m)],color=black):
  od:
  unassign('j'):
  A:=display([seq(p[i],i=0..m)],scaling=constrained,axes=none,thickness=1,tickmarks=[0,0,0]):
  B:=display([seq(q[j],j=0..n)],scaling=constrained,axes=none,thickness=1,tickmarks=[0,0,0]):
  C:=pointplot3d(pts,thickness=4,color=red,tickmarks=[0,0,0]):

  RETURN (display({A,B,C}));
end:

InterPatch:=proc(m,n,Uvals,Vvals,datapts)
  local i,j,k,B,Px,Py,Pz,Bx,By,Bz,U,V,p:
   
  Px:=array(1..m,1..n):
  Py:=array(1..m,1..n):
  Pz:=array(1..m,1..n):
  for i from 1 to m do
    for j from 1 to n do
      k:=i+m*(j-1):
      Px[i,j]:=datapts[k,1]:
      Py[i,j]:=datapts[k,2]:
      Pz[i,j]:=datapts[k,3]:
    od:
  od:
  unassign('i','j','k'):
  
  U:=array(1..m,1..m):
  for i from 1 to m do
    for j from 1 to m do
      B:=binomial(m-1,i-1)*t^(i-1)*(1-t)^(m-i):
      U[i,j]:=subs(t=Uvals[j],B):
    od:
  od: 
  unassign('i','j'): 

  V:=array(1..n,1..n):
  for i from 1 to n do
    for j from 1 to n do
      B:=binomial(n-1,i-1)*t^(i-1)*(1-t)^(n-i):
      V[i,j]:=subs(t=Vvals[j],B):
    od:
  od: 
  unassign('i','j'):

  Bx:=multiply(transpose(inverse(U)),Px,inverse(V)):
  By:=multiply(transpose(inverse(U)),Py,inverse(V));
  Bz:=multiply(transpose(inverse(U)),Pz,inverse(V));
   
  p[1]:=sum(sum(binomial(m-1,i-1)*binomial(n-1,j-1)*u^(i-1)*v^(j-1)*(1-u)^(m-i)*(1-v)^(n-j)*Bx[i,j],i=1..m),j=1..n):
  p[2]:=sum(sum(binomial(m-1,i-1)*binomial(n-1,j-1)*u^(i-1)*v^(j-1)*(1-u)^(m-i)*(1-v)^(n-j)*By[i,j],i=1..m),j=1..n):
  p[3]:=sum(sum(binomial(m-1,i-1)*binomial(n-1,j-1)*u^(i-1)*v^(j-1)*(1-u)^(m-i)*(1-v)^(n-j)*Bz[i,j],i=1..m),j=1..n):

  RETURN([p[1],p[2],p[3]]);
end:

LstSqrFit:=proc(m,n,Uvals,Vvals,datapts,m1,n1)
  local i,j,k,Px,Py,Pz,U,V,UUT,VVT,Bx,By,Bz,p,B:  
   
  Px:=array(1..m,1..n):
  Py:=array(1..m,1..n):
  Pz:=array(1..m,1..n):
  for i from 1 to m do
    for j from 1 to n do
      k:=i+m*(j-1):
      Px[i,j]:=datapts[k,1]:
      Py[i,j]:=datapts[k,2]:
      Pz[i,j]:=datapts[k,3]:
    od:
  od:
  unassign('i','j','k'):
  
  U:=array(1..m1,1..m):
  for i from 1 to m1 do
    for j from 1 to m do
      B:=binomial(m1-1,i-1)*t^(i-1)*(1-t)^(m1-i):
      U[i,j]:=subs(t=Uvals[j],B):
    od:
  od: 
  unassign('i','j'): 

  V:=array(1..n1,1..n):
  for i from 1 to n1 do
    for j from 1 to n do
      B:=binomial(n1-1,i-1)*t^(i-1)*(1-t)^(n1-i):
      V[i,j]:=subs(t=Vvals[j],B):
    od:
  od: 
  unassign('i','j'):
  UUT:=inverse(multiply(U,transpose(U))):
  VVT:=inverse(multiply(V,transpose(V))):
  Bx:=multiply(UUT,U,Px,transpose(V),VVT):
  By:=multiply(UUT,U,Py,transpose(V),VVT);
  Bz:=multiply(UUT,U,Pz,transpose(V),VVT);
   
  p[1]:=sum(sum(binomial(m1-1,i-1)*binomial(n1-1,j-1)*u^(i-1)*v^(j-1)*(1-u)^(m1-i)*(1-v)^(n1-j)*Bx[i,j],i=1..m1),j=1..n1):
  p[2]:=sum(sum(binomial(m1-1,i-1)*binomial(n1-1,j-1)*u^(i-1)*v^(j-1)*(1-u)^(m1-i)*(1-v)^(n1-j)*By[i,j],i=1..m1),j=1..n1):
  p[3]:=sum(sum(binomial(m1-1,i-1)*binomial(n1-1,j-1)*u^(i-1)*v^(j-1)*(1-u)^(m1-i)*(1-v)^(n1-j)*Bz[i,j],i=1..m1),j=1..n1):

  RETURN([p[1],p[2],p[3]]):
end:

SplinePatch:=proc(contnet)
   local numpts,m,n,d,c,b,i,j,J,L,A,B,pts:
   numpts:=rowdim(contnet):
   m:=contnet[numpts,1]:
   n:=contnet[numpts,2]:
   for i from 1 to numpts do
     d[contnet[i,1],contnet[i,2],1]:=contnet[i,3]:
     d[contnet[i,1],contnet[i,2],2]:=contnet[i,4]:
     d[contnet[i,1],contnet[i,2],3]:=contnet[i,5]:
   od:
   for i from 1 to m do
     L:=n-3:     
     c[i,1,1]:=d[i,1,1]: c[i,1,2]:=d[i,1,2]: c[i,1,3]:=d[i,1,3]:
     c[i,2,1]:=d[i,2,1]: c[i,2,2]:=d[i,2,2]: c[i,2,3]:=d[i,2,3]:
     
     c[i,3,1]:=(1/2)*d[i,2,1]+(1/2)*d[i,3,1]:
     c[i,3,2]:=(1/2)*d[i,2,2]+(1/2)*d[i,3,2]:
     c[i,3,3]:=(1/2)*d[i,2,3]+(1/2)*d[i,3,3]:

     for j from 1 to L-2 do
       c[i,3*j+2,1]:=(2/3)*d[i,j+2,1]+(1/3)*d[i,j+3,1]:
       c[i,3*j+2,2]:=(2/3)*d[i,j+2,2]+(1/3)*d[i,j+3,2]:
       c[i,3*j+2,3]:=(2/3)*d[i,j+2,3]+(1/3)*d[i,j+3,3]:
       c[i,3*j+3,1]:=(1/3)*d[i,j+2,1]+(2/3)*d[i,j+3,1]:
       c[i,3*j+3,2]:=(1/3)*d[i,j+2,2]+(2/3)*d[i,j+3,2]:
       c[i,3*j+3,3]:=(1/3)*d[i,j+2,3]+(2/3)*d[i,j+3,3]:
     od:
     
     c[i,3*L-1,1]:=(1/2)*d[i,L+1,1]+(1/2)*d[i,L+2,1]:
     c[i,3*L-1,2]:=(1/2)*d[i,L+1,2]+(1/2)*d[i,L+2,2]:
     c[i,3*L-1,3]:=(1/2)*d[i,L+1,3]+(1/2)*d[i,L+2,3]:
  
     c[i,3*L,1]:=d[i,L+2,1]: c[i,3*L,2]:=d[i,L+2,2]: c[i,3*L,3]:=d[i,L+2,3]:
     c[i,3*L+1,1]:=d[i,L+3,1]: c[i,3*L+1,2]:=d[i,L+3,2]: c[i,3*L+1,3]:=d[i,L+3,3]:

     for j from 1 to L-1 do
       c[i,3*j+1,1]:=(1/2)*c[i,3*j,1]+(1/2)*c[i,3*j+2,1]:
       c[i,3*j+1,2]:=(1/2)*c[i,3*j,2]+(1/2)*c[i,3*j+2,2]:
       c[i,3*j+1,3]:=(1/2)*c[i,3*j,3]+(1/2)*c[i,3*j+2,3]:
     od:
   od:

 for j from 1 to 3*n-8 do
     L:=m-3:
     b[1,j,1]:=c[1,j,1]: b[1,j,2]:=c[1,j,2]: b[1,j,3]:=c[1,j,3]:
     b[2,j,1]:=c[2,j,1]: b[2,j,2]:=c[2,j,2]: b[2,j,3]:=c[2,j,3]:
     
     b[3,j,1]:=(1/2)*c[2,j,1]+(1/2)*c[3,j,1]:
     b[3,j,2]:=(1/2)*c[2,j,2]+(1/2)*c[3,j,2]:
     b[3,j,3]:=(1/2)*c[2,j,3]+(1/2)*c[3,j,3]:

     for i from 1 to L-2 do
       b[3*i+2,j,1]:=(2/3)*c[i+2,j,1]+(1/3)*c[i+3,j,1]:
       b[3*i+2,j,2]:=(2/3)*c[i+2,j,2]+(1/3)*c[i+3,j,2]:
       b[3*i+2,j,3]:=(2/3)*c[i+2,j,3]+(1/3)*c[i+3,j,3]:
       b[3*i+3,j,1]:=(1/3)*c[i+2,j,1]+(2/3)*c[i+3,j,1]:
       b[3*i+3,j,2]:=(1/3)*c[i+2,j,2]+(2/3)*c[i+3,j,2]:
       b[3*i+3,j,3]:=(1/3)*c[i+2,j,3]+(2/3)*c[i+3,j,3]:
     od:

     b[3*L-1,j,1]:=(1/2)*c[L+1,j,1]+(1/2)*c[L+2,j,1]:
     b[3*L-1,j,2]:=(1/2)*c[L+1,j,2]+(1/2)*c[L+2,j,2]:
     b[3*L-1,j,3]:=(1/2)*c[L+1,j,3]+(1/2)*c[L+2,j,3]:
     
     b[3*L,j,1]:=c[L+2,j,1]: b[3*L,j,2]:=c[L+2,j,2]: b[3*L,j,3]:=c[L+2,j,3]:
     b[3*L+1,j,1]:=c[L+3,j,1]: b[3*L+1,j,2]:=c[L+3,j,2]: b[3*L+1,j,3]:=c[L+3,j,3]:

     for i from 1 to L-1 do
       b[3*i+1,j,1]:=(1/2)*b[3*i,j,1]+(1/2)*b[3*i+2,j,1]:
       b[3*i+1,j,2]:=(1/2)*b[3*i,j,2]+(1/2)*b[3*i+2,j,2]:
       b[3*i+1,j,3]:=(1/2)*b[3*i,j,3]+(1/2)*b[3*i+2,j,3]:
     od:
   od:

   for i from 1 to m-3 do
      for j from 1 to n-3 do
         pts:=matrix([
             [b[3*(i-1)+1,3*(j-1)+1,1],b[3*(i-1)+1,3*(j-1)+1,2],b[3*(i-1)+1,3*(j-1)+1,3]],
             [b[3*(i-1)+1,3*(j-1)+2,1],b[3*(i-1)+1,3*(j-1)+2,2],b[3*(i-1)+1,3*(j-1)+2,3]],
             [b[3*(i-1)+1,3*(j-1)+3,1],b[3*(i-1)+1,3*(j-1)+3,2],b[3*(i-1)+1,3*(j-1)+3,3]],
             [b[3*(i-1)+1,3*(j-1)+4,1],b[3*(i-1)+1,3*(j-1)+4,2],b[3*(i-1)+1,3*(j-1)+4,3]],
             [b[3*(i-1)+2,3*(j-1)+1,1],b[3*(i-1)+2,3*(j-1)+1,2],b[3*(i-1)+2,3*(j-1)+1,3]],
             [b[3*(i-1)+2,3*(j-1)+2,1],b[3*(i-1)+2,3*(j-1)+2,2],b[3*(i-1)+2,3*(j-1)+2,3]],
             [b[3*(i-1)+2,3*(j-1)+3,1],b[3*(i-1)+2,3*(j-1)+3,2],b[3*(i-1)+2,3*(j-1)+3,3]],
             [b[3*(i-1)+2,3*(j-1)+4,1],b[3*(i-1)+2,3*(j-1)+4,2],b[3*(i-1)+2,3*(j-1)+4,3]],
             [b[3*(i-1)+3,3*(j-1)+1,1],b[3*(i-1)+3,3*(j-1)+1,2],b[3*(i-1)+3,3*(j-1)+1,3]],
             [b[3*(i-1)+3,3*(j-1)+2,1],b[3*(i-1)+3,3*(j-1)+2,2],b[3*(i-1)+3,3*(j-1)+2,3]],
             [b[3*(i-1)+3,3*(j-1)+3,1],b[3*(i-1)+3,3*(j-1)+3,2],b[3*(i-1)+3,3*(j-1)+3,3]],
             [b[3*(i-1)+3,3*(j-1)+4,1],b[3*(i-1)+3,3*(j-1)+4,2],b[3*(i-1)+3,3*(j-1)+4,3]],
             [b[3*(i-1)+4,3*(j-1)+1,1],b[3*(i-1)+4,3*(j-1)+1,2],b[3*(i-1)+4,3*(j-1)+1,3]],
             [b[3*(i-1)+4,3*(j-1)+2,1],b[3*(i-1)+4,3*(j-1)+2,2],b[3*(i-1)+4,3*(j-1)+2,3]],
             [b[3*(i-1)+4,3*(j-1)+3,1],b[3*(i-1)+4,3*(j-1)+3,2],b[3*(i-1)+4,3*(j-1)+3,3]],
             [b[3*(i-1)+4,3*(j-1)+4,1],b[3*(i-1)+4,3*(j-1)+4,2],b[3*(i-1)+4,3*(j-1)+4,3]]]):
         A[j]:=plot3d(BPatch(3,3,pts),u=0..1,v=0..1,scaling=constrained,style=wireframe):  
      od:
      B[i]:=display([seq(A[j],j=1..n-3)],insequence=false):
   od:

   display({seq(B[i],i=1..m-3)});
end: 


SplineNetA:=proc(contnet)
   local numpts,m,n,d,c,b,L,i,j,A,B,Ac,Ab,pts,P,p,q:
   numpts:=rowdim(contnet):
   m:=contnet[numpts,1]:
   n:=contnet[numpts,2]:
  for i from 1 to numpts do
     d[contnet[i,1],contnet[i,2],1]:=contnet[i,3]:
     d[contnet[i,1],contnet[i,2],2]:=contnet[i,4]:
     d[contnet[i,1],contnet[i,2],3]:=contnet[i,5]:
   od:
   
   for i from 1 to m do
     p[i]:=curve([seq([d[i,j,1],d[i,j,2],d[i,j,3]],j=1..n)],color=red,thickness=2):
   od:
   for j from 1 to n do
     q[j]:=curve([seq([d[i,j,1],d[i,j,2],d[i,j,3]],i=1..m)],color=red,thickness=2):
   od:
   A:=display([seq(p[i],i=1..m)]):
   B:=display([seq(q[j],j=1..n)]):
      display({A,B});
end:

SplineNet:=proc(contnet)
   local numpts,m,n,d,c,b,L,i,j,A,B,Ac,Ab,pts,P,p,q:
   numpts:=rowdim(contnet):
   m:=contnet[numpts,1]:
   n:=contnet[numpts,2]:
  for i from 1 to numpts do
     d[contnet[i,1],contnet[i,2],1]:=contnet[i,3]:
     d[contnet[i,1],contnet[i,2],2]:=contnet[i,4]:
     d[contnet[i,1],contnet[i,2],3]:=contnet[i,5]:
   od:

   for i from 1 to m do
     L:=n-3:     
     c[i,1,1]:=d[i,1,1]: c[i,1,2]:=d[i,1,2]: c[i,1,3]:=d[i,1,3]:
     c[i,2,1]:=d[i,2,1]: c[i,2,2]:=d[i,2,2]: c[i,2,3]:=d[i,2,3]:
     
     c[i,3,1]:=(1/2)*d[i,2,1]+(1/2)*d[i,3,1]:
     c[i,3,2]:=(1/2)*d[i,2,2]+(1/2)*d[i,3,2]:
     c[i,3,3]:=(1/2)*d[i,2,3]+(1/2)*d[i,3,3]:

     for j from 1 to L-2 do
       c[i,3*j+2,1]:=(2/3)*d[i,j+2,1]+(1/3)*d[i,j+3,1]:
       c[i,3*j+2,2]:=(2/3)*d[i,j+2,2]+(1/3)*d[i,j+3,2]:
       c[i,3*j+2,3]:=(2/3)*d[i,j+2,3]+(1/3)*d[i,j+3,3]:
       c[i,3*j+3,1]:=(1/3)*d[i,j+2,1]+(2/3)*d[i,j+3,1]:
       c[i,3*j+3,2]:=(1/3)*d[i,j+2,2]+(2/3)*d[i,j+3,2]:
       c[i,3*j+3,3]:=(1/3)*d[i,j+2,3]+(2/3)*d[i,j+3,3]:
     od:
     
     c[i,3*L-1,1]:=(1/2)*d[i,L+1,1]+(1/2)*d[i,L+2,1]:
     c[i,3*L-1,2]:=(1/2)*d[i,L+1,2]+(1/2)*d[i,L+2,2]:
     c[i,3*L-1,3]:=(1/2)*d[i,L+1,3]+(1/2)*d[i,L+2,3]:
  
     c[i,3*L,1]:=d[i,L+2,1]: c[i,3*L,2]:=d[i,L+2,2]: c[i,3*L,3]:=d[i,L+2,3]:
     c[i,3*L+1,1]:=d[i,L+3,1]: c[i,3*L+1,2]:=d[i,L+3,2]: c[i,3*L+1,3]:=d[i,L+3,3]:

     for j from 1 to L-1 do
       c[i,3*j+1,1]:=(1/2)*c[i,3*j,1]+(1/2)*c[i,3*j+2,1]:
       c[i,3*j+1,2]:=(1/2)*c[i,3*j,2]+(1/2)*c[i,3*j+2,2]:
       c[i,3*j+1,3]:=(1/2)*c[i,3*j,3]+(1/2)*c[i,3*j+2,3]:
     od:
   od:

 for j from 1 to 3*n-8 do
     L:=m-3:
     b[1,j,1]:=c[1,j,1]: b[1,j,2]:=c[1,j,2]: b[1,j,3]:=c[1,j,3]:
     b[2,j,1]:=c[2,j,1]: b[2,j,2]:=c[2,j,2]: b[2,j,3]:=c[2,j,3]:
     
     b[3,j,1]:=(1/2)*c[2,j,1]+(1/2)*c[3,j,1]:
     b[3,j,2]:=(1/2)*c[2,j,2]+(1/2)*c[3,j,2]:
     b[3,j,3]:=(1/2)*c[2,j,3]+(1/2)*c[3,j,3]:

     for i from 1 to L-2 do
       b[3*i+2,j,1]:=(2/3)*c[i+2,j,1]+(1/3)*c[i+3,j,1]:
       b[3*i+2,j,2]:=(2/3)*c[i+2,j,2]+(1/3)*c[i+3,j,2]:
       b[3*i+2,j,3]:=(2/3)*c[i+2,j,3]+(1/3)*c[i+3,j,3]:
       b[3*i+3,j,1]:=(1/3)*c[i+2,j,1]+(2/3)*c[i+3,j,1]:
       b[3*i+3,j,2]:=(1/3)*c[i+2,j,2]+(2/3)*c[i+3,j,2]:
       b[3*i+3,j,3]:=(1/3)*c[i+2,j,3]+(2/3)*c[i+3,j,3]:
     od:

     b[3*L-1,j,1]:=(1/2)*c[L+1,j,1]+(1/2)*c[L+2,j,1]:
     b[3*L-1,j,2]:=(1/2)*c[L+1,j,2]+(1/2)*c[L+2,j,2]:
     b[3*L-1,j,3]:=(1/2)*c[L+1,j,3]+(1/2)*c[L+2,j,3]:
     
     b[3*L,j,1]:=c[L+2,j,1]: b[3*L,j,2]:=c[L+2,j,2]: b[3*L,j,3]:=c[L+2,j,3]:
     b[3*L+1,j,1]:=c[L+3,j,1]: b[3*L+1,j,2]:=c[L+3,j,2]: b[3*L+1,j,3]:=c[L+3,j,3]:

     for i from 1 to L-1 do
       b[3*i+1,j,1]:=(1/2)*b[3*i,j,1]+(1/2)*b[3*i+2,j,1]:
       b[3*i+1,j,2]:=(1/2)*b[3*i,j,2]+(1/2)*b[3*i+2,j,2]:
       b[3*i+1,j,3]:=(1/2)*b[3*i,j,3]+(1/2)*b[3*i+2,j,3]:
     od:
   od:

   for i from 1 to m-3 do
      for j from 1 to n-3 do
         pts:=matrix([
             [b[3*(i-1)+1,3*(j-1)+1,1],b[3*(i-1)+1,3*(j-1)+1,2],b[3*(i-1)+1,3*(j-1)+1,3]],
             [b[3*(i-1)+1,3*(j-1)+2,1],b[3*(i-1)+1,3*(j-1)+2,2],b[3*(i-1)+1,3*(j-1)+2,3]],
             [b[3*(i-1)+1,3*(j-1)+3,1],b[3*(i-1)+1,3*(j-1)+3,2],b[3*(i-1)+1,3*(j-1)+3,3]],
             [b[3*(i-1)+1,3*(j-1)+4,1],b[3*(i-1)+1,3*(j-1)+4,2],b[3*(i-1)+1,3*(j-1)+4,3]],
             [b[3*(i-1)+2,3*(j-1)+1,1],b[3*(i-1)+2,3*(j-1)+1,2],b[3*(i-1)+2,3*(j-1)+1,3]],
             [b[3*(i-1)+2,3*(j-1)+2,1],b[3*(i-1)+2,3*(j-1)+2,2],b[3*(i-1)+2,3*(j-1)+2,3]],
             [b[3*(i-1)+2,3*(j-1)+3,1],b[3*(i-1)+2,3*(j-1)+3,2],b[3*(i-1)+2,3*(j-1)+3,3]],
             [b[3*(i-1)+2,3*(j-1)+4,1],b[3*(i-1)+2,3*(j-1)+4,2],b[3*(i-1)+2,3*(j-1)+4,3]],
             [b[3*(i-1)+3,3*(j-1)+1,1],b[3*(i-1)+3,3*(j-1)+1,2],b[3*(i-1)+3,3*(j-1)+1,3]],
             [b[3*(i-1)+3,3*(j-1)+2,1],b[3*(i-1)+3,3*(j-1)+2,2],b[3*(i-1)+3,3*(j-1)+2,3]],
             [b[3*(i-1)+3,3*(j-1)+3,1],b[3*(i-1)+3,3*(j-1)+3,2],b[3*(i-1)+3,3*(j-1)+3,3]],
             [b[3*(i-1)+3,3*(j-1)+4,1],b[3*(i-1)+3,3*(j-1)+4,2],b[3*(i-1)+3,3*(j-1)+4,3]],
             [b[3*(i-1)+4,3*(j-1)+1,1],b[3*(i-1)+4,3*(j-1)+1,2],b[3*(i-1)+4,3*(j-1)+1,3]],
             [b[3*(i-1)+4,3*(j-1)+2,1],b[3*(i-1)+4,3*(j-1)+2,2],b[3*(i-1)+4,3*(j-1)+2,3]],
             [b[3*(i-1)+4,3*(j-1)+3,1],b[3*(i-1)+4,3*(j-1)+3,2],b[3*(i-1)+4,3*(j-1)+3,3]],
             [b[3*(i-1)+4,3*(j-1)+4,1],b[3*(i-1)+4,3*(j-1)+4,2],b[3*(i-1)+4,3*(j-1)+4,3]]]):
         Ac[j]:=ContNet(3,3,pts):
      od:
      Ab[i]:=display([seq(Ac[j],j=1..n-3)],insequence=false):
   od:
   P:=display([seq(Ab[i],i=1..m-3)],insequence=false):
   display({P});
end: