Generate a Kirkman triple system

Given a universe of $$\v\$$ elements, a Kirkman triple system is a set of $$\(v-1)/2\$$ classes each having $$\v/3\$$ blocks each having three elements, so that

• every pair of elements appears in exactly one block
• all classes are partitions of the universe.

Kirkman's schoolgirl problem corresponds to the $$\v=15\$$ case.

Fifteen young ladies in a school walk out three abreast for seven days in succession: it is required to arrange them daily so that no two shall walk twice abreast.

Below is a procedure to construct a Kirkman triple system for $$\v=3q\$$ where $$\q\$$ is a prime number* of the form $$\6t+1\$$, from my MSE answer here:

Label elements as $$\(x,j)\$$ where $$\x\in\mathbb F_q\$$ and $$\j\in\{0,1,2\}\$$. Let $$\g\$$ be a primitive element of $$\\mathbb F_q\$$. Define blocks $$Z=\{(0,0),(0,1),(0,2)\}\\ B_{i,j}=\{(g^i,j),(g^{i+2t},j),(g^{i+4t},j)\},0\le i and the class $$C=\{Z\}\cup\{B_{i,j}:0\le i Define shifting a block $$\b\$$ by $$\s\in\mathbb F_q\$$ as $$b+s=\{(x+s,j):(x,j)\in b\}$$ and shifting a class similarly, then a Kirkman triple system of order $$\3q\$$ is $$\{C+s:s\in\mathbb F_q\}\cup\{\{A_i+s:s\in\mathbb F_q\}:0\le i<6t,\lfloor i/t\rfloor\in\{0,2,4\}\}$$

Given a prime number $$\q\$$ of the form $$\6t+1\$$, output all classes and blocks of a Kirkman triple system on $$\v=3q\$$ elements. You may use any distinct values for the elements. Formatting is flexible, but the boundaries between elements, blocks and classes must be clear.

This is ; fewest bytes wins. You must be able to run your code to completion for at least the smallest case $$\q=7\$$.

Test cases

This is a possible output for $$\q=7\$$:

[[[0, 7, 14],[1, 2, 4],[8, 9, 11],[15, 16, 18],[3, 13, 19],[6, 12, 17],[5, 10, 20]],
[[1, 8, 15],[2, 3, 5],[9, 10, 12],[16, 17, 19],[4, 7, 20],[0, 13, 18],[6, 11, 14]],
[[2, 9, 16],[3, 4, 6],[10, 11, 13],[17, 18, 20],[5, 8, 14],[1, 7, 19],[0, 12, 15]],
[[3, 10, 17],[0, 4, 5],[7, 11, 12],[14, 18, 19],[6, 9, 15],[2, 8, 20],[1, 13, 16]],
[[4, 11, 18],[1, 5, 6],[8, 12, 13],[15, 19, 20],[0, 10, 16],[3, 9, 14],[2, 7, 17]],
[[5, 12, 19],[0, 2, 6],[7, 9, 13],[14, 16, 20],[1, 11, 17],[4, 10, 15],[3, 8, 18]],
[[6, 13, 20],[0, 1, 3],[7, 8, 10],[14, 15, 17],[2, 12, 18],[5, 11, 16],[4, 9, 19]],
[[1, 9, 18],[2, 10, 19],[3, 11, 20],[4, 12, 14],[5, 13, 15],[6, 7, 16],[0, 8, 17]],
[[2, 11, 15],[3, 12, 16],[4, 13, 17],[5, 7, 18],[6, 8, 19],[0, 9, 20],[1, 10, 14]],
[[4, 8, 16],[5, 9, 17],[6, 10, 18],[0, 11, 19],[1, 12, 20],[2, 13, 14],[3, 7, 15]]]


*The construction also works for $$\q\$$ any prime power of the form $$\6t+1\$$, but I know some languages may be disadvantaged in implementing general finite field arithmetic. Cf. here.

Python with numpy, 285 bytes

from numpy import*
def k(q):
t=q//6;a=array;r=range;R=a(r(q-1));X=R//t%2;h=a(r(3));g=1;p=lambda x:[(x+s)%q for s in r(q)]
while 1in g**R[1:]%q:g+=1
A=a([g**(i+2*h*t)%q+h/5for i in R]);return[p(i)for i in A[X<1]]+p(a([h/5]+[(g**(i+2*h*t))%q+j/5for i in r(t)for j in h]+list(A[X>0])))

Attempt This Online!

PARI/GP, 218 bytes

q->c=concat;c([c([[[[s+0*g=znprimroot(q),j]|j<-r=[0..2]]],c([[[[g^(i+2*t*k)+s,j]|k<-r]|i<-[1..t=q\6]]|j<-r]),[[[g^(i+2*t*j)+s,j]|j<-r]|i<-u=[1..6*t],i\t%2]])|s<-v=[1..q]],[[[[g^(i+2*t*j)+s,j]|j<-r]|s<-v]|i<-u,i\t%2<1])

Attempt This Online!