I tried your method, but it seemed that the expected results could not be produced:

M1={{0,0,1,1/4},{1,0,0,1/4},{0,-1,0,1/4},{0,0,0,1}};
M2={{0,0,-1,0},{0,-1,0,0},{1,0,0,0},{0,0,0,1}};
mod=1
gens={M1,M2};

This is your method:

In[42]:= Union[Join[#1, 
          Flatten[Outer[OperatorApplied[f, {2, 3, 1}][mod], gens, #1], 2]
          (*Flatten[Outer[ AffModOneDotOnLeft, gens, #1, {mod}, 1], 2]*)
          ]] &@{AffineTransform[{IdentityMatrix[3],ConstantArray[0, 3]}]//TransformationMatrix}

Out[42]= {{{{f[-1, 1, 1], f[-1, 0, 1], f[-1, 0, 1], 
    f[-1, 0, 1]}, {f[-1, 0, 1], f[-1, 1, 1], f[-1, 0, 1], 
    f[-1, 0, 1]}, {f[-1, 0, 1], f[-1, 0, 1], f[-1, 1, 1], 
    f[-1, 0, 1]}, {f[-1, 0, 1], f[-1, 0, 1], f[-1, 0, 1], 
    f[-1, 1, 1]}}}, {{{f[0, 1, 1], f[0, 0, 1], f[0, 0, 1], 
    f[0, 0, 1]}, {f[0, 0, 1], f[0, 1, 1], f[0, 0, 1], 
    f[0, 0, 1]}, {f[0, 0, 1], f[0, 0, 1], f[0, 1, 1], 
    f[0, 0, 1]}, {f[0, 0, 1], f[0, 0, 1], f[0, 0, 1], 
    f[0, 1, 1]}}}, {{{f[1/4, 1, 1], f[1/4, 0, 1], f[1/4, 0, 1], 
    f[1/4, 0, 1]}, {f[1/4, 0, 1], f[1/4, 1, 1], f[1/4, 0, 1], 
    f[1/4, 0, 1]}, {f[1/4, 0, 1], f[1/4, 0, 1], f[1/4, 1, 1], 
    f[1/4, 0, 1]}, {f[1/4, 0, 1], f[1/4, 0, 1], f[1/4, 0, 1], 
    f[1/4, 1, 1]}}}, {{{f[1, 1, 1], f[1, 0, 1], f[1, 0, 1], 
    f[1, 0, 1]}, {f[1, 0, 1], f[1, 1, 1], f[1, 0, 1], 
    f[1, 0, 1]}, {f[1, 0, 1], f[1, 0, 1], f[1, 1, 1], 
    f[1, 0, 1]}, {f[1, 0, 1], f[1, 0, 1], f[1, 0, 1], 
    f[1, 1, 1]}}}, {{1, 0, 0, 0}, {0, 1, 0, 0}, {0, 0, 1, 0}, {0, 0, 
   0, 1}}}

This is the desired result:

In[44]:= Union[Join[#1, 
              (*Flatten[Outer[OperatorApplied[f, {2, 3, 1}][mod], 
    gens, #1], 2]*)
              Flatten[Outer[ f, gens, #1, {mod}, 1], 2]
              ]] &@{AffineTransform[{IdentityMatrix[3], 
     ConstantArray[0, 3]}] // TransformationMatrix}

Out[44]= {f[{{0, 0, -1, 0}, {0, -1, 0, 0}, {1, 0, 0, 0}, {0, 0, 0, 
    1}}, {{1, 0, 0, 0}, {0, 1, 0, 0}, {0, 0, 1, 0}, {0, 0, 0, 1}}, 1],
  f[{{0, 0, 1, 1/4}, {1, 0, 0, 1/4}, {0, -1, 0, 1/4}, {0, 0, 0, 
    1}}, {{1, 0, 0, 0}, {0, 1, 0, 0}, {0, 0, 1, 0}, {0, 0, 0, 1}}, 
  1], {{1, 0, 0, 0}, {0, 1, 0, 0}, {0, 0, 1, 0}, {0, 0, 0, 1}}}

Hi Hans,

In fact, I want to keep all the sublists in their original form, i.e., without been flattened. So, I just tweaked your above idea a bit as follows:

In[224]:= 
Join[{x}, Flatten[Outer[f, {a, b}, {x}, {myarg}, 1], 2]]
% /. {x -> Array[x, {2, 2}], a -> Array[a, {1, 2}]}

Out[224]= {x, f[a, x, myarg], f[b, x, myarg]}

Out[225]= {{{x[1, 1], x[1, 2]}, {x[2, 1], x[2, 2]}}, 
 f[{{a[1, 1], a[1, 2]}}, {{x[1, 1], x[1, 2]}, {x[2, 1], x[2, 2]}}, 
  myarg], f[b, {{x[1, 1], x[1, 2]}, {x[2, 1], x[2, 2]}}, myarg]}

Zhao, I have no good answers. In this case I just work from experience and intuition. Guessing, in other words. TreeForm can help in this process. So I leave it to others to explain how Flatten actually handles the level specification {1, 3}.

Check this:

Outer[f, {a, b}, {x, y}, {myarg}] // Flatten[#, {1, 3}] &