Foam Evaluation Program
David Boozer
27 January 2019

Definitions

(* names for vertices of webs *)
v1 = 101;
v2 = 102;
v3 = 103;
v4 = 104;
v5 = 105;
v6 = 106;
v7 = 107;
v8 = 108;
v9 = 109;
v10 = 110;
v11 = 111;
v12 = 112;
v13 = 113;
v14 = 114;
v15 = 115;
v16 = 116;
v17 = 117;
v18 = 118;
v19 = 119;
v20 = 120;
v21 = 121;
v22 = 122;
v23 = 123;
v24 = 124;
v25 = 125;
v26 = 126;
v27 = 127;
v28 = 128;
v29 = 129;
v30 = 130;
v31 = 131;
v32 = 132;
v33 = 133;
v34 = 134;
v35 = 135;

(* names for edges of webs *)
e1 = 201;
e2 = 202;
e3 = 203;
e4 = 204;
e5 = 205;
e6 = 206;
e7 = 207;
e8 = 208;
e9 = 209;
e10 = 210;
e11 = 211;
e12 = 212;
e13 = 213;
e14 = 214;
e15 = 215;
e16 = 216;
e17 = 217;
e18 = 218;
e19 = 219;
e20 = 220;
e21 = 221;
e22 = 222;
e23 = 223;
e24 = 224;
e25 = 225;
e26 = 226;
e27 = 227;
e28 = 228;
e29 = 229;
e30 = 230;
e31 = 231;
e32 = 232;
e33 = 233;
e34 = 234;
e35 = 235;
e36 = 236;
e37 = 237;
e38 = 238;
e39 = 239;
e40 = 240;
e41 = 241;
e42 = 242;
e43 = 243;
e44 = 244;
e45 = 245;
e46 = 246;
e47 = 247;
e48 = 248;
e49 = 249;
e50 = 250;
e51 = 251;
e52 = 252;

(* names for faces of webs *)
f1 = 301;
f2 = 302;
f3 = 303;
f4 = 304;
f5 = 305;
f6 = 306;
f7 = 307;
f8 = 308;
f9 = 309;
f10 = 310;
f11 = 311;
f12 = 312;
f13 = 313;
f14 = 314;
f15 = 315;
f16 = 316;
f17 = 317;
f18 = 318;
f19 = 319;
f20 = 320;

(* names for elementary cobordisms *)
bord1 = 1;
bord2 = 2;
bord3 = 3;
bord4 = 4;
bordZip = 5;
bordIH = 6;
bordSaddle = 7;
bordUnzip = 8;

Define some example webs and foams

web1 :=
{{v1,v2},
{{e1, {v1,v2}},
{e2,{v1,v2}},
{e3,{v1,v2}}},
{{f1,{e1,e2}},
{f2,{e2,e3}}}}

web2 :=
{{v1,v2,v3,v4},
{{e1,{v1,v2}},
{e3,{v1,v2}},
{e4, {v1,v3}},
{e5, {v3,v4}},
{e6,{v3,v4}},
{e7,{v2,v4}}},
{{f1,{e1,e4,e5,e7}},
{f2,{e5,e6}},
{f3,{e3,e4,e6,e7}}}}

web3 :=
{{v1,v2,v3,v4},
{{e1, {v1,v2}},
{e2, {v1,v2}},
{e3, {v1,v4}},
{e4, {v2,v3}},
{e5, {v3,v4}},
{e6, {v3,v4}}},
{{f1, {e1,e2}},
{f2, {e2,e4,e5,e3}},
{f3, {e5,e6}}}}

web4 :=
{{v1, v2, v3, v4, v5, v6},
{{e1, {v1,v2}},
{e2, {v1,v2}},
{e3,{v1,v3}},
{e4, {v2,v4}},
{e5,{v3,v4}},
{e6, {v3,v5}},
{e7, {v4,v6}},
{e8, {v5,v6}},
{e9, {v5,v6}}},
{{f1, {e1,e2}},
{f2, {e2,e4,e5,e3}},
{f3, {e5,e7,e8,e6}},
{f4, {e8,e9}}}}

web5 :=
{{v1, v2, v3, v4, v5, v6, v7, v8},
{{e1, {v1,v2}},
{e2,{v2,v3}},
{e3,{v3,v4}},
{e4,{v4,v1}},
{e5,{v1,v5}},
{e6,{v2,v6}},
{e7,{v3,v7}},
{e8,{v4,v8}},
{e9,{v5,v6}},
{e10,{v6,v7}},
{e11,{v7,v8}},
{e12,{v8,v5}}},
{{f1,{e1,e6,e9,e5}},
{f2,{e2,e7,e10,e6}},
{f3,{e3,e8,e11,e7}},
{f4,{e4,e5,e12,e8}},
{f5,{e9,e10,e11,e12}}}}

web6 := {
{v1,v2,v3,v4},
{{e1,{v1,v2}},
{e2,{v2,v3}},
{e3,{v3,v1}},
{e4,{v1,v4}},
{e5,{v2,v4}},
{e6,{v3,v4}}
},
{{f1,{e1,e5,e4}},
{f2,{e2,e6,e5}},
{f3,{e3,e4,e6}}}}

web7 := {
{v1,v2,v3,v4,v5,v6},
{{e1,{v1,v2}},
{e2,{v2,v3}},
{e3,{v3,v1}},
{e4,{v4,v5}},
{e5,{v5,v6}},
{e6,{v6,v4}},
{e7,{v1,v4}},
{e8,{v2,v5}},
{e9,{v3,v6}}},
{{f1,{e1,e8,e4,e7}},
{f2,{e2,e9,e5,e8}},
{f3,{e3,e7,e6,e9}},
{f4,{e4,e5,e6}}}}

web8 := {
{v1,v2,v3,v4,v5,v6},
{{e1,{v1,v2}},
{e2,{v2,v3}},
{e3,{v3,v4}},
{e4,{v4,v1}},
{e5,{v1,v5}},
{e6,{v2,v6}},
{e7,{v3,v6}},
{e8,{v4,v5}},
{e9,{v5,v6}}},
{{f1,{e1,e6,e9,e5}},
{f2,{e2,e7,e6}},
{f3,{e3,e8,e9,e7}},
{f4,{e4,e5,e8}}}}

webdodec :=
{{v1,v2,v3,v4,v5,v6,v7,v8,v9,v10,v11,v12,v13,v14,v15,v16,v17,v18,v19,v20},
{{e1,{v1,v2}},
{e2,{v2,v3}},
{e3,{v3,v4}},
{e4,{v4,v5}},
{e5,{v5,v1}},
{e6,{v6,v7}},
{e7,{v7,v8}},
{e8,{v8,v9}},
{e9,{v9,v10}},
{e10,{v10,v11}},
{e11,{v11,v12}},
{e12,{v12,v13}},
{e13,{v13,v14}},
{e14,{v14,v15}},
{e15,{v15,v6}},
{e16,{v1,v7}},
{e17,{v2,v9}},
{e18,{v3,v11}},
{e19,{v4,v13}},
{e20,{v5,v15}},
{e21,{v16,v17}},
{e22,{v17,v18}},
{e23,{v18,v19}},
{e24,{v19,v20}},
{e25,{v20,v16}},
{e26,{v6,v16}},
{e27,{v8,v17}},
{e28,{v10,v18}},
{e29,{v12,v19}},
{e30,{v14,v20}}},
{{f1,{e1,e2,e3,e4,e5}},
{f2,{e1,e16,e7,e8,e17}},
{f3,{e2,e17,e9,e10,e18}},
{f4,{e3,e18,e11,e12,e19}},
{f5,{e4,e19,e13,e14,e20}},
{f6,{e5,e20,e15,e6,e16}},
{f7,{e6,e26,e21,e27,e7}},
{f8,{e8,e27,e22,e28,e9}},
{f9,{e10,e28,e23,e29,e11}},
{f10,{e12,e29,e24,e30,e13}},
{f11,{e14,e30,e25,e26,e15}}}
}

foam1 :={
{{},
{}},

{bord1,
{},
{{{e1},2},
{{e2},2},
{{e3},2}}},

{{v1,v2},
{{e1,{v1,v2}, 0},
{e2,{v1,v2}, 0},
{e3,{v1,v2}, 2}}},

{bord2,
{{e3,{e3,e4}, -2}},
{{{e5},2},
{{e6},2}}},

{{v1,v2,v3,v4},
{{e1,{v1,v2}, 0},
{e2,{v1,v2}, 0},
{e3, {v1,v3}, 0},
{e4, {v2,v4}, 2},
{e5,{v3,v4}, 0},
{e6,{v3,v4}, 0}}}
}

NullQ[lst_] :=
SameQ[lst,{}]

NotNullQ[lst_] :=
Not[NullQ[lst]];

IntersectQ[lst1_,lst2_] :=
NotNullQ[Intersection[lst1,lst2]]

SetsEqualQ[lst1_,lst2_] :=
Sort[lst1]==Sort[lst2]

Second[lst_] := lst[[2]]
Third[lst_] := lst[[3]]

(* Example *)
SetsEqualQ[{1,2},{2,1}]
IntersectQ[{1,5},{2,3}]

DisjointQ[lst1_,lst2_] :=
NullQ[Intersection[lst1,lst2]]

SubsetQ[biglst_,smalllst_] :=
SetsEqualQ[Intersection[smalllst,biglst],smalllst]

(* Example *)
DisjointQ[{1,2,3},{3,4}]
DisjointQ[{},{}]
DisjointQ[{},{1}]
SubsetQ[{1,2,3},{}]

PrintLst[lst_] :=
For[i=1,i≤Length[lst],i++,
Print[i," ",lst[[i]]]]

PrintLst[foam1]

TransF[elt_,alist_] :=
Module[
{jlst = Position[Map[First,alist],elt]},
If[NullQ[jlst],elt,Second[alist[[jlst[[1]][[1]]]]]]]

TransF3[elt_,alist_]:=If[alist=={},elt,If[elt==First[First[alist]],First[alist][[3]],TransF3[elt,Rest[alist]]]]

(* Example *)
TransF[1,{{1,11},{2,12},{3,13}}]
TransF3[1,{{1,11,21},{2,12,22},{3,13,23}}]

CheckWeb[web_] :=
Module[
{vertices = First[web],
edges = Map[First,Second[web]],
faces = Map[First,Third[web]]},

Apply[And,
Join[
Map[
Function[x,
Length[Position[
Flatten[Map[Second,Second[web]]],x]]==3],
vertices],
Map[
Function[x,
MemberQ[{1,2},Length[Position[
Flatten[Map[Second,Third[web]]],x]]]],
edges],
{SetsEqualQ[edges,Apply[Union,Map[Second,Third[web]]]]},
{SetsEqualQ[vertices,Apply[Union,Map[Second,Second[web]]]]}
]]]

(* Example *)
CheckWeb[webdodec]

Construct Generating Foams

GetFaceData[face_,web_] :=
Module[
{intvertices,allvertices,intedges,alledges,alledgelst},

intedges = TransF[face,Third[web]];
intvertices =Apply[Union,Map[Function[x,TransF[x,Second[web]]],intedges]];
alledgelst = Select[Second[web],
Function[x,IntersectQ[Second[x],intvertices]]];
alledges = Map[First,alledgelst];
allvertices = Apply[Union,Map[Second,alledgelst]];

{"intvertices"->intvertices,
"extvertices"->Complement[allvertices,intvertices],
"allvertices"-> allvertices,
"intedges"->intedges,
"extedges"->Complement[alledges,intedges],
"alledges"->alledges}
]

(* Example *)
GetFaceData[f1,web3]

DottedWeb[web_,dotlst_] :=
Module[
{dottededges = Map[First,dotlst]},
{
First[web],
Map[Function[x,
If[MemberQ[dottededges,First[x]],
Append[x,TransF[First[x],dotlst]],
Append[x,0]]],
Second[web]]}]

(* Example *)
web1
DottedWeb[web1,{{e1,3}}]

Replace3to1[lst_,old_,new_] :=
Prepend[Drop[RotateLeft[lst,Position[lst,old][[1]][[1]]-2],3],new]

RotFirst[lst_,elt_] :=
RotateLeft[lst,Position[lst,elt][[1]][[1]]-1]

GetPrev[lst_,elt_]:=First[RotateLeft[lst,Position[lst,elt][[1]][[1]]-2]]

GetNext[lst_,elt_]:=First[RotateLeft[lst,Position[lst,elt][[1]][[1]]]]

(* Example *)
Replace3to1[{1, 2, 3, 4, 5,6},1,72]
RotFirst[{1,2,3,4,5},4]
GetNext[{1,2,3,4},2]
GetPrev[{1,2,3,4},1]

VertPrev[edge_,facelst_,web_] :=
Module[
{edgeprev},
If[Length[Second[facelst]]≤2,
Print["VertPrev: not enough edges"];
Abort[]];
edgeprev = Last[RotFirst[Second[facelst],edge]];
First[Intersection[TransF[edge,Second[web]],
TransF[edgeprev,Second[web]]]]]

VertNext[edge_,facelst_,web_] :=
Module[
{edgenext},
If[Length[Second[facelst]]≤2,
Print["VertNext: not enough edges"];
Abort[]];
edgenext =Second[RotFirst[Second[facelst],edge]];
First[Intersection[TransF[edge,Second[web]],
TransF[edgenext,Second[web]]]]]

GetOppFacelst[edge_,facelst_,web_] :=
Module[
{edgeprev,oppedgenext},
edgeprev = Last[RotFirst[Second[facelst],edge]];
oppedgenext = ThirdEdge[edge,edgeprev,web];

temp = Select[Third[web],
Function[x,SubsetQ[Second[x],{edge,oppedgenext}]]];

If[NullQ[temp],{{}},First[temp]]
]

(* Example *)
VertNext[
e9,
First[Select[Third[web5],Function[x,First[x]==f1]]],
web5]
GetOppFacelst[
e9,
First[Select[Third[web5],Function[x,First[x]==f1]]],
web5]

ThirdEdge[edge1_,edge2_,web_] :=
Module[
{vertex = First[Intersection[TransF[edge1,Second[web]],TransF[edge2,Second[web]]]]},
First[First[Select[Second[web],
Function[x,
And[edge1 ≠ First[x], edge2 ≠ First[x],
MemberQ[Second[x],vertex]]]]]]]

(* Example *)
ThirdEdge[e6,e9,web5]
ThirdEdge[e5,e8,web7]

ElimBigon[face_,web_] :=
Module[
{facedata = GetFaceData[face,web],
intvertices,extvertices,allvertices,
intedges,extedges,alledges,
edgeL,edgeR,edgeN,edgeS,edgeO,
bord,newweb},

intvertices = "intvertices"/.facedata;
extvertices = "extvertices"/.facedata;
allvertices = "allvertices"/.facedata;
intedges = "intedges"/.facedata;
extedges = "extedges"/.facedata;
alledges = "alledges"/.facedata;

edgeL = intedges[[1]];
edgeR = intedges[[2]];

If[Length[extedges]==1,

If[Length[DeleteDuplicates[allvertices]] ≠ 2,
Print["ElimBigon: vertices not distinct"];
Abort[]
];

If[Length[DeleteDuplicates[alledges]] ≠ 3,
Print["ElimBigon: edges not distinct"];
Abort[]
];

edgeO  = extedges[[1]];

bord =
{bord1,
{},
{{{edgeR},2},
{{edgeL},2},
{{edgeO},2}}};

newweb = 
{Complement[First[web],intvertices],

Select[Second[web],Function[x,Not[MemberQ[alledges,First[x]]]]],

Select[Third[web],
Function[x,Not[SubsetQ[alledges,Second[x]]]]]};

If[Not[CheckWeb[newweb]],
Print["Invalid web generated by ElimBigon"];
Abort[]];

{{newweb,bord,DottedWeb[web,{{edgeR,0},{edgeL,0}}]},
{newweb,bord,DottedWeb[web,{{edgeR,0},{edgeL,1}}]},
{newweb,bord,DottedWeb[web,{{edgeR,0},{edgeL,2}}]},
{newweb,bord,DottedWeb[web,{{edgeR,1},{edgeL,0}}]},
{newweb,bord,DottedWeb[web,{{edgeR,1},{edgeL,1}}]},
{newweb,bord,DottedWeb[web,{{edgeR,1},{edgeL,2}}]}},

(* Length[extedges] = 2 *)

If[Length[DeleteDuplicates[allvertices]] ≠ 4,
Print["ElimBigon: vertices not distinct"];
Abort[]
];

If[Length[DeleteDuplicates[alledges]] ≠ 4,
Print["ElimBigon: edges not distinct"];
Abort[]
];

edgeN= extedges[[1]];
edgeS = extedges[[2]];

bord =
{bord2,
{{edgeN,{edgeN,edgeS},-2}},
{{{edgeR},2},
{{edgeL},2}}};

newweb = 
{Complement[First[web],intvertices],

Join[{{edgeN,extvertices}},
Select[Second[web],Function[x,Not[MemberQ[alledges,First[x]]]]]],

Map[Function[x,{
First[x],
Which[
MemberQ[Second[x],edgeL],
Replace3to1[Second[x],edgeL,edgeN],
MemberQ[Second[x],edgeR],
Replace3to1[Second[x],edgeR,edgeN],
True,
Second[x]
]}],
Select[Third[web],Function[x,First[x]≠face]]]};

If[Not[CheckWeb[newweb]],
Print["Invalid web generated by ElimBigon"];
Abort[]];

{{newweb,bord,DottedWeb[web,{{edgeL,0}}]},
{newweb,bord,DottedWeb[web,{{edgeL,1}}]}}]
]

(* Example *)
PrintLst[ElimBigon[f1,web1][[5]]]
Print[]
PrintLst[ElimBigon[f2,web2][[1]]]

BordSquare[edgeE_,edgeW_,face_,web_,flagverbose_:False] :=
Module[
{edgeN,edgeS,edgeNE,edgeNW,edgeSW,edgeSE,
newedgeE,newedgeW,
vertNEint,vertNWint,vertSEint,vertSWint,
vertNEext,vertNWext,vertSEext,vertSWext,
facelst,facelstN,facelstS,facelstE,facelstW,
newfacelst,newfacelstE,newfacelstW,
temp,newweb},

facelst =
First[Select[Third[web],Function[x,First[x]==face]]];

If[Length[Second[facelst]] ≠ 4,
Print["BordSquare: face is not a square"];
Abort[]];

If[RotFirst[Second[facelst],edgeE][[3]] ≠ edgeW,
Print["BordSquare: edges are adjacent"];
Abort[]];

edgeN = RotFirst[Second[facelst],edgeE][[2]];
edgeS = RotFirst[Second[facelst],edgeE][[4]];

edgeNE = ThirdEdge[edgeN,edgeE,web];
edgeNW = ThirdEdge[edgeN,edgeW,web];
edgeSW = ThirdEdge[edgeS,edgeW,web];
edgeSE = ThirdEdge[edgeS,edgeE,web];

If[Length[DeleteDuplicates[{edgeNE,edgeNW,edgeSW,edgeSE,edgeN,edgeS,edgeE,edgeW}]] ≠ 8,
Print["BordSquare: edges not distinct"];
Abort[]
];

vertNEint =First[Intersection[TransF[edgeN,Second[web]],TransF[edgeE,Second[web]]]];
vertNWint = First[Intersection[TransF[edgeN,Second[web]],TransF[edgeW,Second[web]]]];
vertSEint = First[Intersection[TransF[edgeS,Second[web]],TransF[edgeE,Second[web]]]];
vertSWint = First[Intersection[TransF[edgeS,Second[web]],TransF[edgeW,Second[web]]]];

vertNEext = First[Complement[TransF[edgeNE,Second[web]],{vertNEint}]];
vertNWext = First[Complement[TransF[edgeNW,Second[web]],{vertNWint}]];
vertSEext = First[Complement[TransF[edgeSE,Second[web]],{vertSEint}]];
vertSWext = First[Complement[TransF[edgeSW,Second[web]],{vertSWint}]];

If[Length[DeleteDuplicates[{vertNEint,vertNWint,vertSEint,vertSWint,vertNEext,vertNWext,vertSEext,vertSWext}]] ≠ 8,
Print["BordSquare: vertices not distinct"];
Abort[]
];

facelstE = GetOppFacelst[edgeE,facelst,web];
facelstW = GetOppFacelst[edgeW,facelst,web];
facelstN = GetOppFacelst[edgeN,facelst,web];
facelstS = GetOppFacelst[edgeS,facelst,web];

newedgeE = edgeE;
newedgeW = edgeW;

newfacelstE =
If[NullQ[First[facelstE]],{{}},
{First[facelstE],
Replace3to1[Second[facelstE],edgeE,newedgeE]}];

newfacelstW =
If[NullQ[First[facelstW]],{{}},
{First[facelstW],
Replace3to1[Second[facelstW],edgeW,newedgeW]}];

newfacelst =
If[Or[NullQ[First[facelstN]],NullQ[First[facelstS]]],{{}},
{First[facelst],
Join[
Drop[RotFirst[Second[facelstN],edgeNW],3],
{newedgeW},
Drop[RotFirst[Second[facelstS],edgeSE],3],
{newedgeE}]}];

If[flagverbose,
Print["facelst: ",facelst];
Print["facelstN: ",facelstN];
Print["facelstS: ",facelstS];
Print["facelstE: ",facelstE];
Print["facelstW: ",facelstW];
Print["newfacelstE: ",newfacelstE];
Print["newfacelstW: ",newfacelstW];
Print["newfacelst: ",newfacelst];

Print["edgeE: ",edgeE];
Print["edgeW: ",edgeW];
Print["edgeN: ",edgeN];
Print["edgeS: ",edgeS];
Print["edgeNE: ",edgeNE];
Print["edgeNW: ",edgeNW];
Print["edgeSE: ",edgeSE];
Print["edgeSW: ",edgeSW];
Print["newedgeE: ",newedgeE];
Print["newedgeW: ",newedgeW];
Print["vertNEint: ",vertNEint];
Print["vertNWint: ",vertNWint];
Print["vertSEint: ",vertSEint];
Print["vertSWint: ",vertSWint];
Print["vertNEext: ",vertNEext];
Print["vertNWext: ",vertNWext];
Print["vertSEext: ",vertSEext];
Print["vertSWext: ",vertSWext];
];

newweb =
{Complement[First[web],{vertNEint,vertNWint,vertSEint,vertSWint}],

Join[Select[Second[web],
Function[x,Not[MemberQ[{edgeN,edgeS,edgeE,edgeW,edgeNE,edgeNW,edgeSE,edgeSW},First[x]]]]],
{{newedgeE, {vertNEext,vertSEext}},
{newedgeW, {vertNWext,vertSWext}}
}],

Join[Select[Third[web],
Function[x,Not[MemberQ[{First[facelst],First[facelstN],First[facelstS],First[facelstE],First[facelstW]},First[x]]]]],
If[NullQ[First[newfacelst]],{},{newfacelst}],
If[NullQ[First[newfacelstE]],{},{newfacelstE}],
If[NullQ[First[newfacelstW]],{},{newfacelstW}]]
};

If[Not[CheckWeb[newweb]],
Print["Invalid web generated by BordSquare"];
Abort[]];

{{newweb,

{bord4,
{{edgeE,{edgeNE,edgeSE},-2},
{edgeW,{edgeNW,edgeSW},-2}},
{{{edgeW},2},
{{edgeE},2},
{{edgeN,edgeS},0}}},

DottedWeb[web,{}]
}}]

(* Example *)
BordSquare[e9,e11,f5,web5][[1]][[1]]
BordSquare[e5,e8,f4,web5][[1]][[1]]
BordSquare[e4,e12,f4,web5][[1]][[1]]

ElimSquare[face_,web_] :=
Module[
{edgelst = TransF[face,Third[web]]},
Join[
BordSquare[edgelst[[1]],edgelst[[3]],face,web],
BordSquare[edgelst[[2]],edgelst[[4]],face,web]]]

ElimTriangle[face_,web_] :=
Module[
{facedata = GetFaceData[face,web],
intvertices,extvertices,allvertices,
intedges,extedges,alledges,
edgeB,edgeR,edgeL,edgeBR,edgeRL,edgeLB,vertex,
newweb},

intvertices = "intvertices"/.facedata;
extvertices = "extvertices"/.facedata;
allvertices = "allvertices"/.facedata;
intedges = "intedges"/.facedata;
extedges = "extedges"/.facedata;
alledges = "alledges"/.facedata;

If[Or[Length[DeleteDuplicates[alledges]] ≠ 6],
Print["ElimTriangle: edges not distinct"];
Print[alledges];
Abort[]
];

edgeB = intedges[[1]];
edgeR = intedges[[2]];
edgeL = intedges[[3]];

edgeBR = ThirdEdge[edgeB,edgeR,web];
edgeRL = ThirdEdge[edgeR,edgeL,web];
edgeLB = ThirdEdge[edgeL,edgeB,web];

vertex =First[Intersection[TransF[edgeL,Second[web]],TransF[edgeR,Second[web]]]];

newweb =
{Append[Complement[First[web],intvertices],vertex],

Join[Select[Second[web],
Function[x,Not[MemberQ[alledges,First[x]]]]],
{{edgeBR,
Append[Intersection[TransF[edgeBR,Second[web]],extvertices],vertex]},
{edgeRL,
Append[Intersection[TransF[edgeRL,Second[web]],extvertices],vertex]},
{edgeLB,
Append[Intersection[TransF[edgeLB,Second[web]],extvertices],vertex]}
}],

Map[Function[x,{
First[x],
Which[
MemberQ[Second[x],edgeR],
Select[Second[x],Function[y,y≠edgeR]],
MemberQ[Second[x],edgeL],
Select[Second[x],Function[y,y≠edgeL]],
MemberQ[Second[x],edgeB],
Select[Second[x],Function[y,y≠edgeB]],
True,
Second[x]
]}],
Select[Third[web],Function[x,First[x] != face]]]};

If[Not[CheckWeb[newweb]],
Print["Invalid web generated by ElimTriangle"];
Print[allvertices];
Print[alledges];
Abort[]];

{{newweb,

{bord3,
{{edgeBR,{edgeBR},-1},
{edgeRL,{edgeRL},-1},
{edgeLB,{edgeLB},-1}},
{{{edgeL},1},
{{edgeB},1},
{{edgeR},1}}},

DottedWeb[web,{}]
}}
]

(* Example *)
ElimTriangle[f4,web7][[1]][[1]]
ElimTriangle[f1,web6][[1]][[1]]

BordUnzip[edgeE_,edgeW_,web_,flagverbose_:False] :=
Module[
{facelst,facelstE,facelstW,facelstN,facelstS,newfacelst,
vertNE,vertNW,vertSE,vertSW,vertN,vertS,
edgeNE,edgeNW,edgeSE,edgeSW,edge,
jlst,newweb},

temp = Select[Third[web],Function[x,
And[MemberQ[Second[x],edgeE],MemberQ[Second[x],edgeW]]]];

If[Length[temp] ≠ 1,
If[flagverbose,
Print["BordUnzip: no unique face containing specfied edges"]];
Return[{}];
];

facelst = First[temp];

vertNE = VertNext[edgeE,facelst,web];
vertSE = VertPrev[edgeE,facelst,web];
vertNW= VertPrev[edgeW,facelst,web];
vertSW = VertNext[edgeW,facelst,web];

If[Length[Union[{vertNE,vertSE,vertNW,vertSW}]]≠4,
If[flagverbose,
Print["BordUnzip: vertices not distinct"]];
Return[{}]
];

temp = Max[Append[First[web],{0}]];
vertN = temp+1;
vertS = temp+2;

temp = Max[Append[Map[First,Third[web]],{0}]];
faceN = temp+1;
faceS = temp+2;

edgeNE = edgeE;
edgeNW = edgeW;

temp = Max[Append[Map[First,Second[web]],{0}]];
edgeSE = temp+1;
edgeSW = temp+2;
edge = temp+3;

temp = RotFirst[Second[facelst],edgeE];
jlst = Position[temp,edgeW];
facelstN =
{faceN,Join[Take[temp,{2,jlst[[1]][[1]]-1}],{edgeNW,edgeNE}]};

temp = RotFirst[Second[facelst],edgeW];
jlst = Position[temp,edgeE];
facelstS =
{faceS,Join[Take[temp,{2,jlst[[1]][[1]]-1}],
{edgeSE, edgeSW}]};

temp = GetOppFacelst[edgeE,facelst,web];
facelstE =
If[NullQ[First[temp]],{{}},
{First[temp],
Join[
Rest[RotFirst[Second[temp],edgeE]],
{edgeNE,edge,edgeSE}]}];

temp = GetOppFacelst[edgeW,facelst,web];
facelstW =
If[NullQ[First[temp]],{{}},
{First[temp],
Join[
Rest[RotFirst[Second[temp],edgeW]],
{edgeSW,edge,edgeNW}]}];

If[flagverbose,
Print["facelst: ",facelst];
Print["facelstN: ",facelstN];
Print["facelstS: ",facelstS];
Print["facelstE: ",facelstE];
Print["facelstW: ",facelstW];
Print["edgeNE: ",edgeNE];
Print["edgeNW: ",edgeNW];
Print["edgeSE: ",edgeSE];
Print["edgeSW: ",edgeSW];
Print["edge: ",edge];
Print["vertNE: ", vertNE];
Print["vertNW: ", vertNW];
Print["vertSE: ", vertSE];
Print["vertSW: ", vertSW];
Print["vertN: ",vertN];
Print["vertS: ",vertS];
];

If[And[NullQ[First[facelstE]],NullQ[First[facelstW]]],
If[flagverbose,
Print["BordUnzip: no E or W face"]];
Return[{}];
];

newweb =
{Join[First[web],{vertN,vertS}],

Join[Select[Second[web],
Function[x,Not[MemberQ[{edgeE,edgeW},First[x]]]]],
{{edgeNE,{vertNE,vertN}}},
{{edgeNW,{vertNW,vertN}}},
{{edgeSE,{vertSE,vertS}}},
{{edgeSW,{vertSW,vertS}}},
{{edge,{vertN,vertS}}}
],

Join[
Select[Third[web],
Function[x,Not[MemberQ[{First[facelst],First[facelstE],First[facelstW]},First[x]]]]],
Select[{facelstE,facelstW,facelstN,facelstS},
Function[x,NotNullQ[First[x]]]]
]};

If[Not[CheckWeb[newweb]],
Print["Invalid web generated by BordUnzip"];
Abort[]];

{{newweb,

{bordUnzip,
{{{edgeNE, edgeSE},{edgeE},-2},
{{edgeNW, edgeSW},{edgeW},-2}},
{{edge,{},2}}},

DottedWeb[web,{}]
}}
]

FaceBordUnzip[face_,web_] :=
Module[
{edgelst = TransF[face,Third[web]],
i,j,foamlst},

foamlst = {};
For[i=1,i≤Length[edgelst],i++,
For[j=i+2,j≤Length[edgelst]-If[i==1,1,0],j++,
foamlst = Join[foamlst,
FullExtendFoamLst[
BordUnzip[edgelst[[i]],edgelst[[j]],web]]];
]
];
foamlst
]

(* Example *)
BordUnzip[e2,e5,web3][[1]][[1]]
BordUnzip[e1,e9,web5][[1]][[1]]

BordSaddle[edgeE_,edgeW_,web_,flagverbose_:False] :=
Module[
{facelst,facelstE,facelstW,facelstN,facelstS,newfacelst,
vertNE,vertNW,vertSE,vertSW,
jlst,newweb},

temp = Select[Third[web],Function[x,
And[MemberQ[Second[x],edgeE],MemberQ[Second[x],edgeW]]]];

If[Length[temp] ≠ 1,
If[flagverbose,
Print["BordSaddle: no unique face containing specified edges"]];
Return[{}]
];

facelst = First[temp];

vertNE = VertNext[edgeE,facelst,web];
vertSE = VertPrev[edgeE,facelst,web];
vertNW= VertPrev[edgeW,facelst,web];
vertSW = VertNext[edgeW,facelst,web];

If[Length[Union[{vertNE,vertSE,vertNW,vertSW}]]≠4,
If[flagverbose,
Print["BordSaddle: vertices not distinct"]];
Return[{}]
];

temp = Max[Append[Map[First,Third[web]],{0}]];
faceN = temp+1;
faceS = temp+2;

temp = RotFirst[Second[facelst],edgeE];
jlst = Position[temp,edgeW];
facelstN =
{faceN,Append[Take[temp,{2,jlst[[1]][[1]]-1}],edgeE]};

temp = RotFirst[Second[facelst],edgeW];
jlst = Position[temp,edgeE];
facelstS =
{faceS,Append[Take[temp,{2,jlst[[1]][[1]]-1}],edgeW]};

facelstE = GetOppFacelst[edgeE,facelst,web];
facelstW = GetOppFacelst[edgeW,facelst,web];

newfacelst =
If[Or[NullQ[First[facelstE]],NullQ[First[facelstW]]],
{{}},
{First[facelst],
Join[
Rest[RotFirst[Second[facelstE],edgeE]],
{edgeE},
Rest[RotFirst[Second[facelstW],edgeW]],
{edgeW}]}];

If[flagverbose,
Print["facelst: ",facelst];
Print["facelstN: ",facelstN];
Print["facelstS: ",facelstS];
Print["facelstE: ",facelstE];
Print["facelstW: ",facelstW];
Print["newfacelst: ",newfacelst];
Print["vertNE: ", vertNE];
Print["vertSE: ", vertSE];
Print["vertNW: ", vertNW];
Print["vertSW: ", vertSW];
];

newweb =
{First[web],

Join[Select[Second[web],
Function[x,Not[MemberQ[{edgeE,edgeW},First[x]]]]],
{{edgeE,{vertNE,vertNW}}},
{{edgeW,{vertSE,vertSW}}}
],

Join[
Select[Third[web],
Function[x,Not[MemberQ[{First[facelst],First[facelstE],First[facelstW]},First[x]]]]],
Select[{newfacelst,facelstN,facelstS},
Function[x,NotNullQ[First[x]]]]
]};

If[Not[CheckWeb[newweb]],
Print["Invalid web generated by BordSaddle"];
Abort[]];

{{newweb,

{bordSaddle,
{{{edgeE, edgeW},{edgeE,edgeW},-4}},
{{}}},

DottedWeb[web,{}]
}}
]

FaceBordSaddle[face_,web_] :=
Module[
{edgelst = TransF[face,Third[web]],
i,j,foamlst},

foamlst = {};
For[i=1,i≤Length[edgelst],i++,
For[j=i+2,j≤Length[edgelst]-If[i==1,1,0],j++,
(*
If[i≠j,Print[edgelst[[i]], ",",edgelst[[j]]]]
*)
foamlst = Join[foamlst,
BordSaddle[edgelst[[i]],edgelst[[j]],web]];
]
];
foamlst
]

(* Example *)
Length[FaceBordSaddle[f5,web5]]
BordSaddle[e2,e5,web3][[1]][[1]]
BordSaddle[e3,e4,web3][[1]][[1]]

BordIH[edge_,web_,flagverbose_:False] :=
Module[
{facelstE,facelstW,facelstN,facelstS,
vertN,vertS,vertE,vertW,
vertNE,vertNW,vertSE,vertSW,
edgeNE,edgeNW,edgeSW,edgeSE,
facelsts,newweb,temp},

temp = Select[Third[web],Function[x,MemberQ[Second[x],edge]]];
facelstE = First[temp];
facelstW = If[Length[temp]==2,Second[temp],{{}}];

edgeNE = GetPrev[Second[facelstE],edge];
edgeSE = GetNext[Second[facelstE],edge];

edgeNW = ThirdEdge[edgeNE,edge,web];
edgeSW = ThirdEdge[edgeSE,edge,web];

vertN =First[Intersection[TransF[edgeNE,Second[web]],TransF[edgeNW,Second[web]]]];

vertS =First[Intersection[TransF[edgeSE,Second[web]],TransF[edgeSW,Second[web]]]];

vertNE = First[Complement[TransF[edgeNE,Second[web]],{vertN}]];
vertNW = First[Complement[TransF[edgeNW,Second[web]],{vertN}]];
vertSE = First[Complement[TransF[edgeSE,Second[web]],{vertS}]];
vertSW = First[Complement[TransF[edgeSW,Second[web]],{vertS}]];

If[Length[Union[{vertNE,vertSE,vertNW,vertSW,vertN,vertS}]]≠6,
If[flagverbose,
Print["BordIH: vertices not distinct"]];
Return[{}]
];

If[Length[Union[{edgeNE,edgeSE,edgeNW,edgeSW,edge}]]≠5,
If[flagverbose,
Print["BordIH: edges not distinct"]];
Return[{}]
];

vertE = vertN;
vertW = vertS;

facelstN = GetOppFacelst[edgeNE,facelstE,web];
facelstS = GetOppFacelst[edgeSE,facelstE,web];

If[And[NullQ[First[facelstN]],NullQ[First[facelstS]]],
If[flagverbose, Print["BordIH: no N or S face"]];
Return[{}]];

If[flagverbose,
Print["facelstN: ",facelstN];
Print["facelstS: ",facelstS];
Print["facelstE: ",facelstE];
Print["facelstW: ",facelstW];
Print["vertN: ",vertN];
Print["vertS: ",vertS];
Print["vertE: ",vertE];
Print["vertW: ",vertW];
Print["vertNE: ",vertNE];
Print["vertNW: ",vertNW];
Print["vertSE: ",vertSE];
Print["vertSW: ",vertSW];
Print["edgeNE: ",edgeNE];
Print["edgeNW: ",edgeNW];
Print["edgeSE: ",edgeSE];
Print["edgeSW: ",edgeSW];
];

newweb=
{First[web],

Join[Select[Second[web],
Function[x,Not[MemberQ[{edge,edgeNE,edgeNW,edgeSE,edgeSW},First[x]]]]],
{{edgeNE,{vertE,vertNE}},
{edgeSE,{vertE,vertSE}},
{edgeNW,{vertW,vertNW}},
{edgeSW,{vertW,vertSW}},
{edge,{vertE,vertW}}
}],

Join[
Select[Third[web],
Function[x,Not[MemberQ[{First[facelstN],First[facelstS],First[facelstE],First[facelstW]},First[x]]]]],
If[NullQ[First[facelstE]],{},
{{First[facelstE],
Drop[RotFirst[Second[facelstE],edge],1]}}],
If[NullQ[First[facelstW]],{},
{{First[facelstW],
Drop[RotFirst[Second[facelstW],edge],1]}}],
If[NullQ[First[facelstN]],{},
{{First[facelstN],
Insert[RotFirst[Second[facelstN],edgeNE],edge,1]}}],
If[NullQ[First[facelstS]],{},
{{First[facelstS],
Insert[RotFirst[Second[facelstS],edgeSW],edge,1]}}]
]};

If[Not[CheckWeb[newweb]],
Print["Invalid web generated by BordIH"];
Abort[]];

{{newweb,

{bordIH,
{{edgeNE,{edgeNE},-1},
{edgeNW,{edgeNW},-1},
{edgeSE,{edgeSE},-1},
{edgeSW,{edgeSW},-1},
{edge, {},1}},
{{{edge},1}}},

DottedWeb[web,{}]
}}
]

(* Example *)
BordIH[e9,web5][[1]][[1]]
BordIH[e3,web3]
BordIH[e5,web4]
BordIH[e6,web5][[1]][[1]]

BordZip[edge_,web_,flagverbose_:False] :=
Module[
{facelstE,facelstW,facelstN,facelstS,
newfacelstE,newfacelstW,facelst,
edgeNE,edgeSE,edgeNW,edgeSW,edgeE,edgeW,
facelsts,temp,newweb,
vertN,vertS,vertNE,vertNW,vertSE,vertSW},

temp = Select[Third[web],Function[x,MemberQ[Second[x],edge]]];

facelstE = First[temp];
facelstW = If[Length[temp]==2,Second[temp],{{}}];

edgeNE = GetPrev[Second[facelstE],edge];
edgeSE = GetNext[Second[facelstE],edge];
edgeNW = ThirdEdge[edgeNE,edge,web];
edgeSW = ThirdEdge[edgeSE,edge,web];

edgeE = edgeNE;
edgeW = edgeNW;

vertN = First[Intersection[TransF[edgeNE,Second[web]],
TransF[edgeNW,Second[web]]]];
vertS = First[Intersection[TransF[edgeSE,Second[web]],
TransF[edgeSW,Second[web]]]];
vertNE = First[Complement[TransF[edgeNE,Second[web]],{vertN}]];
vertNW = First[Complement[TransF[edgeNW,Second[web]],{vertN}]];
vertSE = First[Complement[TransF[edgeSE,Second[web]],{vertS}]];
vertSW =
First[Complement[TransF[edgeSW,Second[web]],{vertS}]];

If[Length[Union[{vertNE,vertSE,vertNW,vertSW,vertN,vertS}]]≠6,
If[flagverbose,
Print["BordZip: vertices not distinct"]];
Return[{}]
];

If[Length[Union[{edgeNE,edgeSE,edgeNW,edgeSW,edge}]]≠5,
If[flagverbose,
Print["BordZip: vertices not distinct"]];
Return[{}]
];

facelstN = GetOppFacelst[edgeNE,facelstE,web];
facelstS = GetOppFacelst[edgeSE,facelstE,web];

newfacelstE =
If[NullQ[First[facelstE]],{{}},
{First[facelstE],Replace3to1[Second[facelstE],edge,edgeNE]}];

If[And[First[facelstN] == First[facelstS],NotNullQ[First[facelstN]]],
If[flagverbose,
Print["BordZip: N and S faces are indentical"]];
Return[{}]];

newfacelstW =
If[NullQ[First[facelstW]],{{}},
{First[facelstW],Replace3to1[Second[facelstW],edge,edgeNW]}];

facelst =
If[Or[NullQ[First[facelstN]],NullQ[First[facelstS]]],
{{}},
{First[facelstN],
Join[{edgeNW},
Drop[RotFirst[Second[facelstS],edgeSE],2],
{edgeNE},
Drop[RotFirst[Second[facelstN],edgeNW],2]]}];

If[flagverbose,
Print["facelstN: ",facelstN];
Print["facelstS: ",facelstS];
Print["facelstE: ",facelstE];
Print["facelstW: ",facelstW];
Print["newfacelstE: ",newfacelstE];
Print["newfacelstW: ",newfacelstW];
Print["facelst: ",facelst];

Print["edgeNE: ",edgeNE];
Print["edgeNW: ",edgeNW];
Print["edgeSE: ",edgeSE];
Print["edgeSW: ",edgeSW];
Print["edgeE: ",edgeE];
Print["edgeW: ",edgeW];

Print["vertN: ",vertN];
Print["vertS: ",vertS];
Print["vertNE: ",vertNE];
Print["vertNW: ",vertNW];
Print["vertSE: ",vertSE];
Print["vertSW: ",vertSW];
];

newweb =
{Complement[First[web],{vertN,vertS}],

Join[Select[Second[web],
Function[x,Not[MemberQ[{edge,edgeNE,edgeNW,edgeSE,edgeSW},First[x]]]]],
{{edgeE,{vertNE,vertSE}},
{edgeW,{vertNW,vertSW}}
}],

Join[
Select[Third[web],
Function[x,Not[MemberQ[{First[facelstN],First[facelstS],First[facelstE],First[facelstW]},First[x]]]]],
If[NullQ[First[facelstE]],{},{newfacelstE}],
If[NullQ[First[facelstW]],{},{newfacelstW}],
If[NullQ[First[facelst]],{},{facelst}]]
};

If[Not[CheckWeb[newweb]],
Print["Invalid web generated by BordZip"];
Abort[]];

{{newweb,

{bordZip,
{{edgeE,{edgeNE,edgeSE},-2},
{edgeW,{edgeNW,edgeSW},-2}},
{{{edge},2}}},

DottedWeb[web,{}]
}}
]

(* Example *)
BordZip[e9,web5][[1]][[1]]
BordZip[e3,web5][[1]][[1]]
BordZip[e8,web5][[1]][[1]]
BordZip[e1,web2]

FindnGon[web_,n_] :=
Module[
{ngonlst = 
Select[Third[web],Function[x,Length[Second[x]]==n]]},
If[NullQ[ngonlst],{},First[First[ngonlst]]]]

(* Example *)
FindnGon[web2,2]
FindnGon[web7,3]

ExtendFoam[foam_] :=
Module[
{web = First[foam],
postfoam = Drop[foam,1],
foamlst},

foamlst= Which[
IntegerQ[FindnGon[web,2]],
ElimBigon[FindnGon[web,2],web],
IntegerQ[FindnGon[web,3]],
ElimTriangle[FindnGon[web,3],web],
IntegerQ[FindnGon[web,4]],
ElimSquare[FindnGon[web,4],web],
True,
(*
Print["No bigons, triangles, or squares"];
*)
Return[{}]];

Map[Function[x,Join[x,postfoam]],foamlst]]

ExtendFoamLst[foamlst_] :=
Apply[Join,
Map[
Function[x,If[NullQ[First[First[x]]],{x},ExtendFoam[x]]],
foamlst]]

FullExtendFoamLst[foamlst_] :=
If[Apply[And,Map[Function[x,NullQ[First[First[x]]]],foamlst]],
foamlst,
FullExtendFoamLst[ExtendFoamLst[foamlst]]]

(* Example *)
Length[FullExtendFoamLst[{{web1}}]]

GetBordFoamlstU[web_]:=Module[{faces=Map[First,Third[web]]},Apply[Join,Map[Function[x,FullExtendFoamLst[FaceBordUnzip[x,web]]],faces]]]

GetBordFoamlstZ[web_]:=Module[{edges=Map[First,Second[web]]},Apply[Join,Map[Function[x,FullExtendFoamLst[BordZip[x,web]]],edges]]]

GetBordFoamlstS[web_]:=Module[{faces=Map[First,Third[web]]},Apply[Join,Map[Function[x,FullExtendFoamLst[FaceBordSaddle[x,web]]],faces]]]

GetBordFoamlstI[web_]:=Module[{edges=Map[First,Second[web]]},Apply[Join,Map[Function[x,FullExtendFoamLst[BordIH[x,web]]],edges]]]

GetBordFoamlst[web_] :=
Module[
{edges = Map[First,Second[web]],
faces = Map[First,Third[web]],
lstUnzip,lstZip,lstSaddle,lstIH,n},

lstZip =
Apply[Join,Map[
Function[x,FullExtendFoamLst[BordZip[x,web]]],
edges]];

lstIH =
Apply[Join,Map[
Function[x,FullExtendFoamLst[BordIH[x,web]]],
edges]];

lstSaddle =
Apply[Join,Map[
Function[x,FullExtendFoamLst[FaceBordSaddle[x,web]]],faces]];

lstUnzip =
Apply[Join,Map[
Function[x,FullExtendFoamLst[FaceBordUnzip[x,web]]],
faces]];

n = Length[lstUnzip] + Length[lstZip] + Length[lstSaddle] + Length[lstIH];

Print["U: ", Length[lstUnzip], 
" (", Length[lstUnzip]/n//N, ")"];
Print["Z: ", Length[lstZip], 
" (", Length[lstZip]/n//N, ")"];
Print["I: ", Length[lstIH], 
" (", Length[lstIH]/n//N, ")"];
Print["S: ", Length[lstSaddle], 
" (", Length[lstSaddle]/n//N, ")"];
Join[lstUnzip,lstZip,lstSaddle,lstIH]
];

Bordnum2Type[n_]:=TransF[n,{{bordIH,"I"},{bordSaddle,"S"},{bordUnzip,"U"},{bordZip,"Z"}}]

Foam Evaluator

FlipFoam[foam_] :=
Module[
{n = Length[foam]},
Table[foam[[n-i+1]],{i,1,n}]]

JoinFacets[facetlst_] :=
If[NullQ[facetlst],{},
Module[
{totallst = Transpose[facetlst]},
Table[
If[OddQ[i],
Apply[Join,totallst[[i]]],
{0,Apply[Plus,Map[Second,totallst[[i]]]]}
],{i,1,Length[totallst]}]]]

Boundary[edgelst_,web_] :=
DeleteDuplicates[Flatten[Map[Function[x,TransF[x,Second[web]]],edgelst]]]

(* Example *)
Boundary[{e6,e9},web5]
facet1 = {{},{1,2},{201},{0,0},{201}};
facet2 = {{},{1,2},{202},{0,0},{202}};
JoinFacets[{facet1,facet2}]

FacetlstDisjointQ[facetlst_,foam_] :=
Module[
{n,temp,i},
temp = DeleteDuplicates[Map[Length,facetlst]];
If[Length[temp]≠1,
Print["FacetlstDisjointQ: facets have unequal lengths"];
Abort[];
];
n = First[temp];
For[i=1,i≤n,i+=2,
temp = 
Map[Function[x,Boundary[x[[i]],foam[[i]]]],facetlst];
(*
Print["i=",i," temp=",temp];
*)
If[Length[Apply[Union,temp]]<Length[Apply[Join,temp]],
Return[False]]];
True]

(* example *)
facet1 = {{},{0,0},{e1},{0,0},{e1}};facet2 = {{},{0,0},{},{0,0},{e5}};
FacetlstDisjointQ[{facet1,facet2},foam1]

facet1 = {{},{0,0},{e1}};
facet2 = {{},{0,0},{e3}};
FacetlstDisjointQ[{facet1,facet2},foam1]

ReplaceEdges[edgelst_,bordlst_,bordtype_] :=
Module[
{jlst,i,flag,newedgelst,addlst,euler},

flag = False;
newedgelst = edgelst;
euler = 0;
addlst = {};
For[i=1,i≤Length[bordlst],i++,
jlst = Position[newedgelst,First[bordlst[[i]]]];
If[NotNullQ[jlst],
flag = True;
newedgelst = Delete[newedgelst,jlst[[1]][[1]]];
addlst = Join[addlst,Second[bordlst[[i]]]];
euler = euler + Third[bordlst[[i]]];
]
];
If[flag,
{{bordtype,euler},Join[newedgelst,addlst]},
{{0,0},edgelst}]
]

(* Example *)
bordlst = {{203,{203,204},-2},{205,{207,208},-2}};
ReplaceEdges[{e1,e2,e3,e5},bordlst,3]

CheckPartition[P_,S_] :=
Module[
{union = Apply[Union,P]},
And[SetsEqualQ[union,S],
Apply[Plus,Map[Length,P]]==Length[S]]]

(* Example *)
CheckPartition[{{1,2},{2,3}},{1,2,3}]
CheckPartition[{{1,2},{3}},{1,2,3}]

CreateNullFacet[n_] :=
Table[If[OddQ[i],{},{0,0}],{i,1,n}]

ExtendFacetlst[facetlst_,foam_] :=
Module[
{bord,bordtype,bordlstreplace,bordlstnew,n,
edgeE,edgeW,jE,jW,facetE,facetW,i,
edgeNE,edgeNW,edgeSE,edgeSW,edge,
jNE,jNW,jSE,jSW,j,
facetNE,facetSE,facetNW,facetSW,facet,
flag,flagE,flagW},

n = If[NullQ[facetlst],1,Length[First[facetlst]]];
bord = foam[[n+1]];
bordtype = First[bord];
bordlstreplace = Second[bord];
bordlstnew = Third[bord];

If[Not[CheckPartition[Map[Last,facetlst],Map[First,Second[foam[[n]]]]]],
Print["ExtendFacetLst: invalid facetlst"];
Print[facetlst];
Abort[];
];

Which[
(* Saddle *)
bordtype == bordSaddle,

(* bottom edges *)
edgeE = bordlstreplace[[1]][[1]][[1]];
edgeW = bordlstreplace[[1]][[1]][[2]];
jE = 0;
jW = 0;
For[i=1,i<=Length[facetlst],i++,
If[MemberQ[Last[facetlst[[i]]],edgeE],jE=i];
If[MemberQ[Last[facetlst[[i]]],edgeW],jW=i]
];
facetE = facetlst[[jE]];
facetW = facetlst[[jW]];
flag = Not[FacetlstDisjointQ[{facetE,facetW},foam]];
If[And[flag,jE ≠ jW],
Print["BordSaddle: self-intersection"];
Return[{}];
];
Join[
Map[
Function[x,Join[x,{{0,0}, Last[x]}]],
Delete[facetlst,{{jE},{jW}}]
],
If[flag,
{Join[
facetE,
{{bordSaddle,-4}},
{Last[facetE]}]},
{Join[
JoinFacets[{facetE,facetW}],
{{bordSaddle,-4}},
{Join[Last[facetE],Last[facetW]]}]}
]
],

(* Unzip *)
bordtype == bordUnzip,

(* bottom edges *)
edgeNE = bordlstreplace[[1]][[1]][[1]];
edgeSE = bordlstreplace[[1]][[1]][[2]];
edgeNW = bordlstreplace[[2]][[1]][[1]];
edgeSW = bordlstreplace[[2]][[1]][[2]];
edge = bordlstnew[[1]][[1]];
jNE = 0;
jSE = 0;
jNW = 0;
jSW = 0;
j = 0;
For[i=1,i<=Length[facetlst],i++,
If[MemberQ[Last[facetlst[[i]]],edgeNE],jNE=i];
If[MemberQ[Last[facetlst[[i]]],edgeSE],jSE=i];
If[MemberQ[Last[facetlst[[i]]],edgeNW],jNW=i];
If[MemberQ[Last[facetlst[[i]]],edgeSW],jSW=i];
If[MemberQ[Last[facetlst[[i]]],edge],j=i];
];
facetNE = facetlst[[jNE]];
facetSE = facetlst[[jSE]];
facetNW = facetlst[[jNW]];
facetSW = facetlst[[jSW]];
facet = facetlst[[j]];
flagE = Not[FacetlstDisjointQ[{facetNE,facetSE},foam]];
flagW = Not[FacetlstDisjointQ[{facetNW,facetSW},foam]];
If[Or[And[flagE,jNE ≠ jSE],And[flagW, jNW  ≠ jSW]],
(*
Print["BordUnzip: facet self-intersection"];
*)
Return[{}]
];
Join[
Map[
Function[x,Join[x,{{0,0}, Last[x]}]],
Delete[facetlst,{{jNE},{jSE},{jNW},{jSW},{j}}]
],
If[flagE,
{Join[
facetNE,
{{bordUnzip,-2}},{Append[Complement[Last[facetNE],{edgeNE,edgeSE}],bordlstreplace[[1]][[2]][[1]]]}]},
{Join[
JoinFacets[{facetNE,facetSE}],
{{bordUnzip,-2}},
{Append[
Complement[Join[Last[facetNE],Last[facetSE]],
{edgeNE,edgeSE}],
bordlstreplace[[1]][[2]][[1]]]}]}
],
If[flagW,
{Join[
facetNW,
{{bordUnzip,-2}},
{Append[
Complement[Last[facetNW],{edgeNW,edgeSW}],
bordlstreplace[[2]][[2]][[1]]]}]},
{Join[
JoinFacets[{facetNW,facetSW}],
{{bordUnzip,-2}},
{Append[
Complement[Join[Last[facetNW],Last[facetSW]],
{edgeNW,edgeSW}],
bordlstreplace[[2]][[2]][[1]]]}]}
],
{Join[
facet,
{{bordUnzip, 2}},
{Complement[Last[facet],{edge}]}]}
],

True,
(* any bordtype other than bordSaddle or bordUnzip *)

If[bordtype==bordIH,
(* bottom edges *)
edgeNE = bordlstreplace[[1]][[1]];
edgeNW = bordlstreplace[[2]][[1]];
edgeSE = bordlstreplace[[3]][[1]];
edgeSW = bordlstreplace[[4]][[1]];
If[
Apply[Or,
Map[
Function[x,
Or[
And[MemberQ[Last[x],edgeNE],MemberQ[Last[x],edgeNW]],
And[MemberQ[Last[x],edgeSE],MemberQ[Last[x],edgeSW]]
]],
facetlst]],
(*
Print["BordIH: facet self-intersection"];
*)
Return[{}]
]
];
If[bordtype==bordZip,
(* bottom edges *)
edgeE = bordlstreplace[[1]][[1]];
edgeW = bordlstreplace[[2]][[1]];
If[
Apply[Or,
Map[
Function[x,
And[MemberQ[Last[x],edgeE],MemberQ[Last[x],edgeW]]],
facetlst]],
(*
Print["BordZip: facet self-intersection"];
*)
Return[{}]
]
];
Join[
(* create new facets *)
Map[Function[x,
Join[
CreateNullFacet[n],
{{bordtype,Second[x]}},
{First[x]}]],
bordlstnew],

(* extend existing facets *)
Map[
Function[x,
Join[x,ReplaceEdges[Last[x],bordlstreplace,bordtype]]],
facetlst],
{}]
]]

(* Example *)
PrintLst[ExtendFacetlst[{},foam1]]
Print[]
PrintLst[ExtendFacetlst[ExtendFacetlst[{},foam1],foam1]]

GetFacetsHelper[facetlst_,foam_] :=
Which[
NullQ[facetlst],
{},
Length[First[facetlst]]==Length[foam],
facetlst,
True,
GetFacetsHelper[ExtendFacetlst[facetlst,foam],foam]]

GetFacets[foam_] :=
GetFacetsHelper[ExtendFacetlst[{},foam],foam]

(* Example *)
PrintLst[GetFacets[foam1]]

MergePartitions[P1_,P2_] :=
If[
NullQ[P2],
P1,
MergePartitions[
Append[
Select[P1,Function[x,DisjointQ[x,First[P2]]]],Apply[Union,
Select[P1,Function[x,IntersectQ[x,First[P2]]]]]],
Rest[P2]]]

(* Example *)
MergePartitions[{{1,2},{3,4}},{{1,4},{2,3}}]

FlipFacet[facet_] :=
Module[
{n = Length[facet]},
Table[facet[[n-i+1]],{i,1,n}]]

GlueFacets[facet1_,facet2_] :=
Join[facet1,Drop[FlipFacet[facet2],1]]

PairFacets[foam1_,foam2_,verboseflag_:False] :=
Module[
{facets1,facets2,nullfacets1,nullfacets2,
P1,P2,P,joinfacets1,joinfacets2,pairedfacets,i},

facets1 = GetFacets[foam1];
facets2 = GetFacets[foam2];
If[Or[NullQ[facets1],NullQ[facets2]],
Return[{}]];

nullfacets1 = Select[facets1,Function[x,NullQ[Last[x]]]];
nullfacets2 = Select[facets2,Function[x,NullQ[Last[x]]]];

P1 = Select[Map[Last,facets1],NotNullQ];
P2 = Select[Map[Last,facets2],NotNullQ];
P = MergePartitions[P1,P2];

If[Not[And[
CheckPartition[P1,Map[First,Second[Last[foam1]]]],CheckPartition[P2,Map[First,Second[Last[foam2]]]]]],
Print["PairFacets: invalid facetlst"];
Abort[];
];

If[verboseflag,
Print["P1: ",P1];
Print["S1: ",Map[First,Second[Last[foam1]]]];
Print["P2: ",P2];
Print["S2: ",Map[First,Second[Last[foam2]]]];
Print["P: ",P];
Print["# of nullfacets1: ",Length[nullfacets1]];
Print["# of nullfacets2: ",Length[nullfacets2]];
];

pairedfacets = {};
For[i=1,i≤Length[P],i++,
joinfacets1=
Select[facets1,Function[y,IntersectQ[Last[y],P[[i]]]]];
joinfacets2=
Select[facets2,Function[y,IntersectQ[Last[y],P[[i]]]]];

If[And[
FacetlstDisjointQ[joinfacets1,foam1],
FacetlstDisjointQ[joinfacets2,foam2]],
pairedfacets = 
Append[pairedfacets,
GlueFacets[
JoinFacets[joinfacets1],JoinFacets[joinfacets2]]],
Return[{}]]];

Join[
pairedfacets,
Map[
Function[x,GlueFacets[x,CreateNullFacet[Length[foam2]]]],
nullfacets1],
Map[
Function[x,GlueFacets[CreateNullFacet[Length[foam1]],x]],
nullfacets2]]
]

(* Example of foam with no admissible colorings *)
foamlstdodec = FullExtendFoamLst[{BordZip[e1,webdodec][[1]]}];
testfoam1= foamlstdodec[[13]];
testfoam2 = foamlstdodec[[26]];
Length[testfoam1]
Length[testfoam2]
Length[PairFacets[testfoam1,testfoam2]]

(* Example *)
foam8 =FullExtendFoamLst[BordIH[e9,web8]][[1]];
PrintLst[PairFacets[foam8,foam8]]

(* Example (should be 5) *)
Length[PairFacets[
FullExtendFoamLst[{ElimSquare[f5,web5][[1]]}][[1]],FullExtendFoamLst[{ElimSquare[f5,web5][[2]]}][[1]]]]

(* Example (should be 8) *)
Length[PairFacets[
FullExtendFoamLst[{ElimSquare[f5,web5][[1]]}][[1]],FullExtendFoamLst[{ElimSquare[f5,web5][[1]]}][[1]]]]

(* Example (should be 5) *)
Length[PairFacets[
FullExtendFoamLst[{{web3}}][[1]],FullExtendFoamLst[{{web3}}][[1]]]]

(* Example (should be 6) *)
Length[PairFacets[
FullExtendFoamLst[{{web6}}][[1]],FullExtendFoamLst[{{web6}}][[1]]]]

(* Example (should be 5) *)
Length[PairFacets[
FullExtendFoamLst[{{web2}}][[1]],FullExtendFoamLst[{{web2}}][[1]]]]

(* Example (should be 6) *)
Length[PairFacets[
FullExtendFoamLst[{BordZip[e9,web8][[1]]}][[1]],FullExtendFoamLst[{BordZip[e9,web8][[1]]}][[1]]]]

(* Example (should be 6) *)
Length[PairFacets[
FullExtendFoamLst[{BordZip[e9,web8][[1]]}][[1]],FullExtendFoamLst[{ElimTriangle[f4,web8][[1]]}][[1]]]]

(* glue two half foams to create a closed foam *)
GlueFoams[foam1_,foam2_] :=
Module[
{last1 = Last[foam1],
last2 = Last[foam2]},
Join[
Drop[foam1,-1],
{{First[last1],
Table[
{First[Second[last1][[i]]],
Second[Second[last1][[i]]],
Third[Second[last1][[i]]] + Third[Second[last2][[i]]]},
{i,1,Length[Second[last1]]}]}},
Drop[FlipFoam[foam2],1]]]

MakeTable[facetlst_,foam_]:=Module[{vlist,newfacetlst,table,i,j,k,temp},vlist=Table[If[OddQ[i],First[foam[[i]]],{}],{i,1,Length[foam]}];
newfacetlst=Table[Table[If[OddQ[i],Apply[Plus,Map[Function[x,2^(First[First[Position[vlist[[i]],x]]]-1)],Boundary[facetlst[[j]][[i]],foam[[i]]]]],0],{i,1,Length[foam]}],{j,1,Length[facetlst]}];
table=ConstantArray[True,{Length[facetlst],Length[facetlst]}];
For[i=1,i≤Length[facetlst],i++,For[j=i+1,j≤Length[facetlst],j++,temp=Apply[Or,Table[BitAnd[newfacetlst[[i]][[k]],newfacetlst[[j]][[k]]]≠0,{k,1,Length[foam],2}]];
table[[i]][[j]]=temp;
table[[j]][[i]]=temp;]];
table]

MakeGraph[facetlst_,foam_] :=
Select[Flatten[Table[
If[
Not[FacetlstDisjointQ[{facetlst[[i]],facetlst[[j]]},foam]],
{i,j},
{}],
{i,1,Length[facetlst]},{j,i+1,Length[facetlst]}],1],
Function[x,NotNullQ[x]]]

(* Example *)
foam2 = FullExtendFoamLst[{{web2}}][[1]];
testfoam = GlueFoams[foam2,foam2];
facets2 = PairFacets[foam2,foam2];
MakeTable[facets2,testfoam]//MatrixForm
MakeGraph[facets2,testfoam]

Colorings[vlist_,graph_,disallowedlst_]  :=
If[NullQ[vlist],{{}},
Module[
{v = First[vlist]},
Apply[Join,
Table[
If[
MemberQ[disallowedlst,{v,i}],{},
Map[
Function[x,Join[{{v,i}},x]],
Colorings[
Rest[vlist],
graph,
Join[
Select[
Map[
Function[x,
Which[
First[x]==v,{Second[x],i},
Second[x]==v,{First[x],i},
True,{}]],
graph],
NotNullQ],
disallowedlst]]]],
{i,1,3}]]]]

(* Example *)
PrintLst[Colorings[{1,2,3},{{1,2},{1,3},{2,3}},{}]]

FindAColoring[vlist_,graph_,disallowedlst_,c_]  :=
If[NullQ[vlist],c,
Module[
{v = First[vlist],cnew,i},
For[i=1,i≤3,i++,
If[
Not[MemberQ[disallowedlst,{v,i}]],

cnew =
FindAColoring[
Rest[vlist],
graph,
Join[
Select[
Map[
Function[x,
Which[
First[x]==v,{Second[x],i},
Second[x]==v,{First[x],i},
True,{}]],
graph],
NotNullQ],
disallowedlst],
Append[c,{v,i}]];

If[Not[NullQ[cnew]],Return[cnew]];
]];
{}]]

FindAColoring[{1,2,3},{{1,2},{1,3},{2,3}},{},{}]

SetBit[x_,i_] :=
Switch[i,
1, BitOr[x,1],
2, BitOr[x,2],
3, BitOr[x,4]]

BitSetQ[x_,i_] :=
Switch[i,
1, Not[BitAnd[x,1]==0],
2, Not[BitAnd[x,2]==0],
3, Not[BitAnd[x,4]==0]]

(* Example *)
SetBit[SetBit[0,2],1]
BitSetQ[6,2]

SetElt[lst_,n_,i_] :=
Module[
{newlst = lst},
newlst[[n]] = i;
newlst]

Lst2alst[lst_] :=
Table[{i,lst[[i]]},{i,1,Length[lst]}]

Graph2Table[graph_] :=
Module[
{n = Max[Flatten[graph]],
table,i,j},

table = Table[Table[False,{i,1,n}],{j,1,n}];
For[i=1,i≤Length[graph],i++,
table[[First[graph[[i]]]]][[Second[graph[[i]]]]]=True;
table[[Second[graph[[i]]]]][[First[graph[[i]]]]]=True;
];
table
]

FastColoringsHelper[n_,table_,disallowedlst_]  :=
Block[
{newdislst=disallowedlst,pair,i,k,colorings},

If[n==0,Return[{Table[0,{i,1,Length[disallowedlst]}]}]];
colorings = {};
Apply[Join,
Table[
If[
BitSetQ[disallowedlst[[n]],i],
{},
newdislst=disallowedlst;
For[k=1,k≤Length[table],k++,
If[table[[n]][[k]],
newdislst[[k]] = SetBit[newdislst[[k]],i]]
];

Map[
Function[x,SetElt[x,n,i]],
FastColoringsHelper[n-1,table,newdislst]]],
{i,1,3}]]]

FastColorings[n_,table_,disallowedlst_] :=
Map[Lst2alst,FastColoringsHelper[n,table,disallowedlst]]

Tait[web_] :=
Module[
{vlist,graph},
vlist = Map[First,Second[web]];

graph = Select[Flatten[Table[
If[
IntersectQ[TransF[vlist[[i]],Second[web]],
TransF[vlist[[j]],Second[web]]],
{i,j},
{}],
{i,1,Length[vlist]},{j,i+1,Length[vlist]}],1],
Function[x,NotNullQ[x]]];

Length[FastColorings[Length[vlist],Graph2Table[graph],ConstantArray[0,Length[vlist]]]]
]

(* Example *)
Tait[web5]

Euler[facet_] :=
Sum[Second[facet[[i]]],{i,2,Length[facet]-1,2}]/4

EulerS[c1_,c2_,coloring_,facets_] :=
Module[
{euler=
Sum[
If[Or[TransF[i,coloring]==c1,TransF[i,coloring]==c2],
Euler[facets[[i]]],0],
{i,1,Length[facets]}]},
If[EvenQ[euler],euler,
Print["Invalid Euler characteristic: ",euler];
Abort[]]]

PolyP[facets_,dots_,coloring_]:=
Product[
Switch[
TransF[i,coloring],
1,x1^dots[[i]],
2,x2^dots[[i]],
3,x3^dots[[i]]
],
{i,1,Length[facets]}]

PolyQ[facets_,coloring_]:=(x1+x2)^(EulerS[1,2,coloring,facets]/2)*(x2+x3)^(EulerS[2,3,coloring,facets]/2)*(x3+x1)^(EulerS[3,1,coloring,facets]/2)

CountDots[facet_,foam_]:=
Sum[
Sum[
TransF3[facet[[i]][[j]],Second[foam[[i]]]],{j,1,Length[facet[[i]]]}],
{i,1,Length[foam],2}]

GetDotLst[facets_,foam_]:=Map[Function[x,CountDots[x,foam]],facets]

Evaluate foams in various rings

(* evaluate foam in the ring R = F2[X1,X2,X3] *)
EvalR[foam1_,foam2_]:=
Module[
{facets,table,dotlst,colorings,foam,temp,
t1,t2,t3,t4,t},

{t1,facets}=Timing[PairFacets[foam1,foam2]];
If[NullQ[facets],Return[0]];

foam=GlueFoams[foam1,foam2];
{t2,table}=Timing[MakeTable[facets,foam]];
dotlst=GetDotLst[facets,foam];

{t3,colorings}=
Timing[FastColorings[Length[facets],table,ConstantArray[0,Length[facets]]]];

{t4,eval} = Timing[Together[Sum[PolyP[facets,dotlst,colorings[[i]]]/PolyQ[facets,colorings[[i]]],{i,1,Length[colorings]}]]];

(*
Print["(",t1,", ",t2, ", ", t3, ", ", t4, ")"];
t = t1+t2+t3+t4;
Print["(",t1/t,", ",t2/t, ", ", t3/t, ", ", t4/t, ")"];
*)

temp = PolynomialReduce[Numerator[eval],Denominator[eval],{x1,x2,x3},Modulus->2];
If[temp[[2]]≠0,
Print["EvalR: remainder not zero"];
Abort[];
];
temp[[1]][[1]]
]

(* evaluate foam in the ring F2[E1,E2,E3] *)
EvalSymm[foam1_,foam2_] :=
Module[
{temp =
PolynomialMod[
SymmetricReduction[ EvalR[foam1,foam2],{x1,x2,x3},{s1,s2,s3}],2]},

If[temp[[2]]≠0,
Print["EvalSymm: remainder not zero"];
Abort[];
];
temp[[1]]
]

(* evaluate foam in the ring F2[E] *)
EvalPID[foam1_,foam2_] :=
EvalSymm[foam1,foam2]/.{s1->0,s2->0,s3->x}

(* evaluate foam in the field F2 *)
EvalBit[foam1_,foam2_] :=
EvalR[foam1,foam2]/.{x1->0,x2->0,x3->0}

(* Example *)
foamlstdodec = FullExtendFoamLst[{BordZip[e1,webdodec][[1]]}];
Timing[EvalR[foamlstdodec[[5]],foamlstdodec[[12]]]]
Timing[EvalSymm[foamlstdodec[[5]],foamlstdodec[[12]]]]
Timing[EvalPID[foamlstdodec[[5]],foamlstdodec[[12]]]]
Timing[EvalBit[foamlstdodec[[5]],foamlstdodec[[12]]]]

(* Example *)
Timing[EvalR[foamlstdodec[[9]],foamlstdodec[[12]]]]
Timing[EvalR[foamlstdodec[[6]],foamlstdodec[[9]]]]

Timing[EvalSymm[foamlstdodec[[5]],foamlstdodec[[12]]]]
Timing[EvalSymm[foamlstdodec[[9]],foamlstdodec[[12]]]]
Timing[EvalSymm[foamlstdodec[[6]],foamlstdodec[[9]]]]

GetMatrix[foamlst1_,foamlst2_] :=
Table[EvalBit[foamlst1[[i]],foamlst2[[j]]],
{i,1,Length[foamlst1]},{j,1,Length[foamlst2]}]

(* Example *)
foamlst1 = FullExtendFoamLst[{{web1}}];
Timing[A = GetMatrix[foamlst1,foamlst1];]
A //MatrixForm
MatrixRank[A,Modulus->2]

NumDots[foam_] :=
Apply[Plus,Flatten[Table[Map[Third,Second[foam[[i]]]],
{i,1,Length[foam],2}]]]

Bord2Deg[bordtype_] :=
Which[
bordtype==bord1,
-3,
bordtype == bord2,
-1,
bordtype == bord3,
0,
bordtype==bord4,
0,
bordtype==bordZip,
1,
bordtype==bordUnzip,
1,
bordtype==bordSaddle,
2,
bordtype==bordIH,
1,
True,
Print["Bord2Deg: Error"];
Abort[];
]

DegUndotted[foam_] :=
Apply[Plus,Table[Bord2Deg[First[foam[[i]]]],{i,2,Length[foam],2}]]

Deg[foam_] :=
2*NumDots[foam] + DegUndotted[foam]

foamlst = FullExtendFoamLst[{{web6}}];
U = foamlst[[1]];
Deg[U]

A =Table[EvalSymm[foamlst[[i]],foamlst[[j]]],{i,1,6},{j,1,6}];
A//MatrixForm

Table[Deg[foamlst[[i]]]+Deg[foamlst[[j]]],{i,1,6},{j,1,6}]//MatrixForm

(* evaluate foam in the ring R = F2[X1,X2,X3] *)
FastEvalBit[foam1_,foam2_]:=
Module[
{facets,table,dotlst,colorings,eval,foam,temp,
deg1,deg2,
t1,t2,t3,t4,t},

deg1 = Deg[foam1];
deg2 = Deg[foam2];
If[Or[NullQ[deg1],NullQ[deg2]],Return[0]];
If[deg1 + deg2 ≠ 0,Return[ 0]];

{t1,facets}=Timing[PairFacets[foam1,foam2]];
If[NullQ[facets],Return[0]];

foam=GlueFoams[foam1,foam2];
{t2,table}=Timing[MakeTable[facets,foam]];
dotlst=GetDotLst[facets,foam];

{t3,colorings}=
Timing[FastColorings[Length[facets],table,ConstantArray[0,Length[facets]]]];

{t4,eval} = Timing[Together[Sum[PolyP[facets,dotlst,colorings[[i]]]/PolyQ[facets,colorings[[i]]],{i,1,Length[colorings]}]]];

(*
Print["(",t1,", ",t2, ", ", t3, ", ", t4, ")"];
t = t1+t2+t3+t4;
Print["(",t1/t,", ",t2/t, ", ", t3/t, ", ", t4/t, ")"];
*)
temp = PolynomialReduce[Numerator[eval],Denominator[eval],{x1,x2,x3},Modulus->2];
If[temp[[2]]≠0,
Print["FastEvalBit: remainder not zero"];
Abort[];
];
temp[[1]][[1]]/.{x1->0,x2->0,x3->0}
]

(* Example *)
foamlst1 = FullExtendFoamLst[{{web1}}];
FastEvalBit[foamlst1[[1]],foamlst1[[6]]]
EvalBit[foamlst1[[1]],foamlst1[[6]]]

Linear Algebra

FindBasis[vlist_] :=
Module[
{numvectors = Length[vlist],
A = Transpose[vlist],
B,dim,genlst,i,j},

dim = Length[A];
B = RowReduce[A,Modulus->2];
genlst = {};
For[i=1,i<=dim,i++,
For[j=1,And[j<=numvectors,B[[i]][[j]]==0],j++];
(*
Print["i=", i, " j=",j];
*)
If[j≤numvectors,genlst=Append[genlst,j]];
];
genlst
]

(* Example *)
FindBasis[
{{1,0},
{1,0},
{1,1},
{0,1}}]

PolynomialLCM2[polylst__] :=
PolynomialLCM[polylst,Modulus->2]

DropFirstCol[A_] :=
Module[
{B= Drop[Transpose[A],1]},
If[NullQ[B],{},Transpose[B]]]

(* check if a polynomial f in F2[x] is zero *)
PolyZeroQ[f_] :=
SameQ[PolynomialMod[f,2],0]

MatrixRankPID[A_] :=
Module[
{B=A,f,i,j,k,rank,v},

rank=0;
For[k=1,k≤Length[A[[1]]],k++,
(*
Print[k];
Print[B//MatrixForm];
*)
v=Transpose[B][[1]];
For[i=1,And[i≤Length[v],PolyZeroQ[v[[i]]]],i++];
If[i<=Length[v],
rank++;
f = Apply[PolynomialLCM2,
Select[v,Function[x,Not[PolyZeroQ[x]]]]];
B = Drop[
Table[
If[
PolyZeroQ[v[[j]]],
B[[j]],
PolynomialMod[B[[j]]*PolynomialQuotient[f,v[[j]],x,Modulus->2]-
B[[i]]*PolynomialQuotient[f,v[[i]],x,Modulus->2],2]
],
{j,1,Length[v]}],
{i}]
];
If[Not[NullQ[B]],B = DropFirstCol[B]];
];
rank]

MatrixRankPID[
{{1,0,1},
{1,x,0},
{0,x,1}}]

Smith Normal Form

MatrixMod2[A_]:=
Table[
PolynomialMod[A[[i]][[j]],2],
{i,1,Length[A]},{j,1,Length[A[[1]]]}]

PolyZeroQ[f_]:=
SameQ[PolynomialMod[f,2],0]

PolyZeroQ[2+2 t]

CombineRows[A_,i_,j_,a_,b_,c_,d_]:=
Table[
Which[
k==i,a*A[[i]]+b*A[[j]],
k==j,c*A[[i]]+d*A[[j]],
True,A[[k]]],
{k,1,Length[A]}]

FixRow[A_,S_,i_,j_,col_]:=
Module[
{a=A[[i]][[col]],
b=A[[j]][[col]],d,x,y,ap,bp},

If[PolyZeroQ[b],Return[{A,S}]];

{d,{x,y}}=
If[PolyZeroQ[a-b],
{a,{1,0}},
PolynomialExtendedGCD[a,b,t,Modulus->2]];

(*
Print["FixRow: ",
x*PolynomialQuotient[a,d,t,Modulus->2] +
y*PolynomialQuotient[b,d,t,Modulus->2]];
*)
ap = PolynomialQuotient[a,d,t,Modulus->2];
bp = PolynomialQuotient[b,d,t,Modulus->2];
Simplify[{
CombineRows[A,i,j,x,y,-bp,ap],CombineRows[S,i,j,x,y,-bp,ap]}]]

FixCol[A_,B_,i_,j_,row_]:=
Module[
{Ap,Bp},{Ap,Bp}=FixRow[Transpose[A],Transpose[B],i,j,row];
{Transpose[Ap],Transpose[Bp]}]

SwapRows[A_,i_,j_]:=
Table[
Which[
k==i,A[[j]],
k==j,A[[i]],
True,A[[k]]],
{k,1,Length[A]}]

SwapCols[A_,i_,j_]:=
Transpose[SwapRows[Transpose[A],i,j]];

A={{6,2,3},{4,5,6},{8,8,9}};
A//MatrixForm
SwapCols[A,2,3]//MatrixForm

SmithForm[A_]:=
Module[
{n=Length[A],
m=Length[A[[1]]],
B=A/.x->t,
S,T,i,j,pivotold,k,ip,jp,flag},

S=IdentityMatrix[n];
T=IdentityMatrix[m];

For[k=1,k≤n-1,k++,
flag=True;
For[i=k,And[i≤n,flag],i++,For[j=k,And[j≤m,flag],j++,If[Not[PolyZeroQ[B[[i]][[j]]]],{ip,jp}={i,j};
flag=False]]];

If[flag,Break[]];

B=SwapRows[B,k,ip];
S=SwapRows[S,k,ip];
B=SwapCols[B,k,jp];
T=SwapCols[T,k,jp];

While[True,
pivotold=B[[k]][[k]];
For[i=k+1,i≤n,i++,{B,S}=FixRow[B,S,k,i,k]];
For[i=k+1,i≤m,i++,{B,T}=FixCol[B,T,k,i,k]];
If[PolyZeroQ[B[[k]][[k]]-pivotold],Break[]];
];
];
{MatrixMod2[S],MatrixMod2[B],MatrixMod2[T]}
]

MatrixRankPIDF2[A_]:=
Module[
{B=SmithForm[A/.x->t][[2]],
rankF2=0,rankPID=0,k},

For[k=1,k≤Length[B],k++,If[Not[SameQ[B[[k]][[k]]/.t->0,0]],rankF2++];
If[Not[PolyZeroQ[B[[k]][[k]]]],rankPID++];
];
{rankPID,rankF2}]

A:={{1+x^2,x},{x+1,x^2}}
SmithForm[A][[2]]//MatrixForm

A:={{1+t^2,t},{t+1,t^2}}
SmithForm[A][[2]]//MatrixForm
MatrixRankPIDF2[A]
S := SmithForm[A][[1]];
B := SmithForm[A][[2]];
T := SmithForm[A][[3]];
MatrixMod2[(1+2t)*Inverse[S].B.Inverse[T]]//MatrixForm

A :=(({{0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0}, {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, x}, {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, {1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1}, {0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0}, {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, {0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0}, {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, {0, x, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0}}))
SmithForm[A][[2]]//MatrixForm

Find generators

FindGeneratorsPID[foamlst_,flagparallel_:False] :=
Module[
{A,t1,ttot,i,j,k,v,w,z,b,rank,oldrank,foundlst,genlst,temp,r1,r2,S,B,T,deg,newdeg,degindex},

A = {{EvalPID[foamlst[[1]],foamlst[[1]]]}};

oldrank = MatrixRankPID[A];
foundlst = {{1,oldrank}};
ttot = 0;
For[i=2,i≤Length[foamlst],i++,

If[flagparallel,
w = Take[foamlst,i-1];
z = foamlst[[i]];
{t1,v} = 
Timing[{ParallelTable[EvalPID[z,w[[j]]],
{j,1,i-1},Method->"CoarsestGrained"]}],

{t1,v} = 
Timing[
{Table[EvalPID[foamlst[[i]],foamlst[[j]]],{j,1,i-1}]}]
];

ttot += t1;

b = EvalPID[foamlst[[i]],foamlst[[i]]];

A = ArrayFlatten[{{A,Transpose[v]},{v,b}}];
rank = MatrixRankPID[A];
If[rank > oldrank,
oldrank = rank;
foundlst = Append[foundlst,{i,rank}];
];

Print["i=",i, " rankPID=", rank,
" rankF2=", MatrixRank[A/.x->0,Modulus->2],
" type=",
Bordnum2Type[First[foamlst[[i]][[Length[foamlst[[i]]]-1]]]], " t=",t1, " ttot=", ttot];

If[i==Length[foamlst],

{S,B,T} = SmithForm[A];
r1=0;
r2=0;
For[k=1,k≤Length[B],k++,
If[Not[PolyZeroQ[B[[k]][[k]]]],r1++];
If[Not[SameQ[B[[k]][[k]]/.t->0,0]],r2++];
];
Print["rankPID=", r1, " rankF2=",r2];
temp = {};
For[k=1,k≤r1,k++,
deg = {};
For[j=1,j<=Length[T],j++,
If[Not[PolyZeroQ[T[[j]][[k]]]],
newdeg=
Deg[foamlst[[j]]]+6*Exponent[T[[j]][[k]],t];

If[NullQ[deg],
degindex=j;
deg=newdeg,
If[deg≠newdeg,
Print["Generator not homogeneous: k=",k,
" j=",j, " deg=",deg, " newdeg=",newdeg,
" oldj=",degindex];
Print[T//MatrixForm];
Abort[];
]
]
];
];
temp=Append[temp,deg];
];
If[Not[NullQ[temp]],
Print[Apply[Plus,
Table[q[k]*Count[temp,k],{k,Min[temp],Max[temp]}]]]];

genlst = FindBasis[A/.x->0];
temp =
Table[Deg[foamlst[[genlst[[k]]]]],{k,1,Length[genlst]}];
If[Not[NullQ[temp]],
Print[Apply[Plus,
Table[q[k]*Count[temp,k],{k,Min[temp],Max[temp]}]]]];
]
];
{}]

foamlst1 = FullExtendFoamLst[{{web1}}];
A =FindGeneratorsPID[foamlst1]

FindGenerators[foamlst_,flagparallel_:False] :=
Module[
{A,B,t,ttot,i,j,v,w,z,b,rank,oldrank,foundlst,genlst,temp},

A = {{FastEvalBit[foamlst[[1]],foamlst[[1]]]}};

oldrank = MatrixRank[A,Modulus->2];
foundlst = {{1,oldrank}};
ttot = 0;
For[i=2,i≤Length[foamlst],i++,

If[flagparallel,
w = Take[foamlst,i-1];
z = foamlst[[i]];
{t,v} = 
Timing[{ParallelTable[
FastEvalBit[z,w[[j]]],
{j,1,i-1},Method->"CoarsestGrained"]}],

{t,v} = 
Timing[{Table[
FastEvalBit[foamlst[[i]],foamlst[[j]]],
{j,1,i-1}]}]
];

ttot += t;
b = FastEvalBit[foamlst[[i]],foamlst[[i]]];

A = ArrayFlatten[{{A,Transpose[v]},{v,b}}];
rank = MatrixRank[A,Modulus->2];
If[rank > oldrank,
oldrank = rank;
foundlst = Append[foundlst,{i,rank}];
];
Print["i=",i, " j=", rank, " type=",
Bordnum2Type[First[foamlst[[i]][[Length[foamlst[[i]]]-1]]]], " t=",t, " ttot=", ttot];


genlst = FindBasis[A];
temp =
Table[Deg[foamlst[[genlst[[k]]]]],{k,1,Length[genlst]}];
If[Not[NullQ[temp]],
Print[Apply[Plus,
Table[q[k]*Count[temp,k],{k,Min[temp],Max[temp]}]]]];

If[Mod[i,100]==0,Print[foundlst]];
];
Print[foundlst];
{}]

foamlst1 = FullExtendFoamLst[{{web1}}];
FindGenerators[foamlst1]

SaveGenerators[filename_,foamlst_,flagparallel_:False] :=
Module[
{B,t,ttot,i,j,v,w,z,b,rank,oldrank,foundlst,genlst,temp},

If[Length[A]==0,
Print["Initializing A"];
A = {{FastEvalBit[foamlst[[1]],foamlst[[1]]]}};
];

oldrank = MatrixRank[A,Modulus->2];
foundlst = {{1,oldrank}};
ttot = 0;
For[i=Length[A]+1,i≤Length[foamlst],i++,

If[flagparallel,
w = Take[foamlst,i-1];
z = foamlst[[i]];
{t,v} = 
Timing[{ParallelTable[
FastEvalBit[z,w[[j]]],
{j,1,i-1},Method->"CoarsestGrained"]}],

{t,v} = 
Timing[{Table[
FastEvalBit[foamlst[[i]],foamlst[[j]]],
{j,1,i-1}]}]
];

ttot += t;
b = FastEvalBit[foamlst[[i]],foamlst[[i]]];

A = ArrayFlatten[{{A,Transpose[v]},{v,b}}];
rank = MatrixRank[A,Modulus->2];
If[rank > oldrank,
oldrank = rank;
foundlst = Append[foundlst,{i,rank}];
];
Print["i=",i, " j=", rank, " type=",
Bordnum2Type[First[foamlst[[i]][[Length[foamlst[[i]]]-1]]]], " t=",t, " ttot=", ttot];


genlst = FindBasis[A];
temp =
Table[Deg[foamlst[[genlst[[k]]]]],{k,1,Length[genlst]}];
If[Not[NullQ[temp]],
Print[Apply[Plus,
Table[q[k]*Count[temp,k],{k,Min[temp],Max[temp]}]]]];

If[Mod[i,100]==0,
DeleteFile[filename];
Print["Saving A to ", filename];
Save[filename,A];
Print[foundlst]];
];
Print[foundlst];
{}]