(* Malgorzata Mikosz, Andrzej Weber
Triality in so(4,4). *)

(* Two generic elements of so(4,4) *)
DD=Table[x[i+10j],{j,1,8},{i,1,8}];
Do[DD[[9-i,i]]=0,{i,1,8}];Do[Do[DD[[i,j]]=-DD[[9-j,9-i]],{i,10-j,8}],{j,2,8}]
CC=Table[y[i+10j],{j,1,8},{i,1,8}];
Do[CC[[9-i,i]]=0,{i,1,8}];
Do[Do[CC[[i,j]]=-CC[[9-j,9-i]],{i,10-j,8}],{j,2,8}]

(* The triality automorphism *)
triality[A_]:={{1/2 (A[[1,1]]+A[[2,2]]+A[[3,3]]+A[[4,4]]),A[[3,5]],-A[[2,5]],A[[1,5]],A[[2,6]],A[[1,6]],A[[1,7]],0},{A[[5,3]],1/2 (A[[1,1]]+A[[2,2]]-A[[3,3]]-A[[4,4]]),A[[2,3]],-A[[1,3]],A[[2,4]],A[[1,4]],0,-A[[1,7]]},{-A[[5,2]],A[[3,2]],1/2 (A[[1,1]]-A[[2,2]]+A[[3,3]]-A[[4,4]]),A[[1,2]],A[[3,4]],0,-A[[1,4]],-A[[1,6]]},{A[[5,1]],-A[[3,1]],A[[2,1]],1/2 (-A[[1,1]]+A[[2,2]]+A[[3,3]]-A[[4,4]]),0,-A[[3,4]],-A[[2,4]],-A[[2,6]]},{A[[6,2]],A[[4,2]],A[[4,3]],0,1/2 (A[[1,1]]-A[[2,2]]-A[[3,3]]+A[[4,4]]),-A[[1,2]],A[[1,3]],-A[[1,5]]},{A[[6,1]],A[[4,1]],0,-A[[4,3]],-A[[2,1]],1/2 (-A[[1,1]]+A[[2,2]]-A[[3,3]]+A[[4,4]]),-A[[2,3]],A[[2,5]]},{A[[7,1]],0,-A[[4,1]],-A[[4,2]],A[[3,1]],-A[[3,2]],1/2 (-A[[1,1]]-A[[2,2]]+A[[3,3]]+A[[4,4]]),-A[[3,5]]},{0,-A[[7,1]],-A[[6,1]],-A[[6,2]],-A[[5,1]],A[[5,2]],-A[[5,3]],1/2 (-A[[1,1]]-A[[2,2]]-A[[3,3]]-A[[4,4]])}};
 MatrixForm[triality[DD]]

(* 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 bracket *)
Simplify[bracket[triality[CC],triality[DD]]]==Simplify[triality[bracket[CC,DD]]]

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

HH=GG/.{x[61]->a11,x[71]->a12,x[22]->t2,x[32]->a9,x[52]->a8,x[62]->a10,x[23]->a7,x[33]->t3,x[53]->a6,x[25]->a2,x[35]->a1,x[16]->a4,x[26]->a3,x[17]->a5};MatrixForm[HH]

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

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

(* The basis of (e5-e4)^\perp *)
f[i_]:=(vector=IdentityMatrix[8][[i]];If[i==4,vector[[5]]=1];vector)
(* The vector f[4] is equal to e4+e5 *)

(* The basis of 3 forms *)
base=Complement[Range[1,8],{5}];
base3=Subsets[base,{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[{3,4,6}]]