(* Malgorzata Mikosz, Andrzej Weber
Triality so(8) following E. Cartan *)

(* General matrices from so(8) *)
CC=Table[x[i+10j],{j,1,8},{i,1,8}];
Do[Do[CC[[i,j]]=-CC[[j,i]],{j,1,i}];CC[[i,i]]=0,{i,1,8}]
DD=Table[y[i+10j],{j,1,8},{i,1,8}];
Do[Do[DD[[i,j]]=-DD[[j,i]],{j,1,i}];DD[[i,i]]=0,{i,1,8}];MatrixForm[DD]

(* Rotation matrix *)
Mrot={{-(1/2),-(1/2),-(1/2),-(1/2)},{(1/2),(1/2),-(1/2),-(1/2)},{(1/2),-(1/2),(1/2),-(1/2)},{(1/2),-(1/2),-(1/2),(1/2)}};
rot=Function[{a,b,c,d},Simplify[Mrot.{a,b,c,d}]];

(* The distinguisched quadruple *)
mod[a_]:=Mod[a-1,7]+2
Do[orb[i]={{1,i+1},{mod[i+1],mod[i+5]},{mod[i+4],mod[i+6]},{mod[i+2],mod[i+3]}},{i,1,7}]

(* The definition of triality automorphism *)
triality[a_]:=(bm=Table[0,{i,1,8},{j,1,8}];
Do[Do[bm[[orb[i][[k]][[1]],orb[i][[k]][[2]]]]=rot[a[[orb[i][[1]][[1]],orb[i][[1]][[2]]]],a[[orb[i][[2]][[1]],orb[i][[2]][[2]]]],a[[orb[i][[3]][[1]],orb[i][[3]][[2]]]],a[[orb[i][[4]][[1]],orb[i][[4]][[2]]]]][[k]],{k,1,4}],{i,1,7}];bm=bm-Transpose[bm]
;bm)

(* The order of the triality is equal to 3 *)
Simplify[triality[triality[triality[DD]]]]==DD

(* Definition of the Lie bracket *)
bracket[A_,B_]:=Simplify[A.B-B.A]

(* Triality preserves the bracket *)
bracket[triality[CC],triality[DD]]==Simplify[triality[bracket[CC,DD]]]

(* G2 - The form of matrices fixed by the triality *)
GG=CC/.(SolveAlways[triality[CC]==CC,{}][[1]]);
MatrixForm[GG]

xx=Variables[GG];aa={a1,a2,a3,a4,a5,a6,a7,a8,a9,a10,a11,a12,a13,a14};
HH=GG/.Table[xx[[i]]->aa[[i]],{i,1,14}];MatrixForm[HH]

v={b1,b2,b3,b4,b5,b6,b7,b8};

(* The vector fixed by G2 *)
v/.SolveAlways[{v}.HH==0,Variables[HH]]

(* The basis of (e1)^\perp *)
f[i_]:=IdentityMatrix[8][[i]]

(* The basis of 3 forms *)
base3=Subsets[Range[2,8],{3}]

(* Coefficients in the space of 3-forms *)
wyb[a_,b_]:={b[[a[[1]]]],b[[a[[2]]]],b[[a[[3]]]]}
minor[a_,v1_,v2_,v3_]:=Det[wyb[a,Transpose[{v1,v2,v3}]]]

(* The action of the general matrix *)
h[v1_]:=HH.v1 ;
(* The action on \Lambda^3 *)
h3[v1_,v2_,v3_]:=Factor[Table[minor[a,h[v1],v2,v3]+ minor[a,v1,h[v2],v3]+ minor[a,v1,v2,h[v3]],{a,base3}]]

(* Fixed 3-form *)
variables=Table[c[a],{a,base3}];
virt=Sum[c[a]e[a],{a,base3}];
Factor[virt/.(SolveAlways[Sum[c[a]h3[f[a[[1]]],f[a[[2]]],f[a[[3]]]],{a,base3}]==0,Variables[HH]][[1]])]

Factor[%/ c[{5,7,8}]]

(* Note, that here the basis of im(O) is e2,e3,...,e8. It is shifted by one comparing with Cartan's notation *)