```% Data

% Blocks - point at origin is implicitly included
block([p(1,0,0),p(0,1,0),p(0,0,1)]).
block([p(0,1,0),p(1,1,0),p(1,2,0)]).
block([p(0,1,0),p(1,1,0),p(1,1,1)]).
block([p(0,1,0),p(0,1,1),p(1,1,1)]).
block([p(1,0,0),p(1,1,0),p(2,0,0)]).
block([p(0,1,0),p(1,0,0),p(2,0,0)]).
block([p(1,0,0),p(0,1,0)]).

% 24 posssible permutation of the axes
direction(v(1,2,3)).
direction(v(1,-2,-3)).
direction(v(1,3,-2)).
direction(v(1,-3,2)).
direction(v(-1,2,-3)).
direction(v(-1,-2,3)).
direction(v(-1,3,2)).
direction(v(-1,-3,-2)).
direction(v(2,1,-3)).
direction(v(2,-1,3)).
direction(v(2,3,1)).
direction(v(2,-3,-1)).
direction(v(-2,1,3)).
direction(v(-2,-1,-3)).
direction(v(-2,3,-1)).
direction(v(-2,-3,1)).
direction(v(3,1,2)).
direction(v(3,-1,-2)).
direction(v(3,2,-1)).
direction(v(3,-2,1)).
direction(v(-3,1,-2)).
direction(v(-3,-1,2)).
direction(v(-3,2,1)).
direction(v(-3,-2,-1)).

% Generating rotations

select_axis(1, p(X,_,_), X).
select_axis(-1, p(X,_,_), X1) :- X1 is -X.
select_axis(2, p(_,Y,_), Y).
select_axis(-2, p(_,Y,_), Y1) :- Y1 is -Y.
select_axis(3, p(_,_,Z), Z).
select_axis(-3, p(_,_,Z), Z1) :- Z1 is -Z.

permute_axes(v(I,J,K), Point, p(X,Y,Z)) :-
select_axis(I, Point, X),
select_axis(J, Point, Y),
select_axis(K, Point, Z).

rotation(Block, Rotated) :-
direction(D),
maplist(permute_axes(D), [p(0,0,0)|Block], Rotated).

minimum([Point], Point).
minimum([p(X,Y,Z)|Points], p(X1,Y1,Z1)) :-
minimum(Points, p(X0,Y0,Z0)),
X1 is min(X, X0),
Y1 is min(Y, Y0),
Z1 is min(Z, Z0).

subtract(p(X1,Y1,Z1), p(X2,Y2,Z2), p(X,Y,Z)) :-
X is X2 - X1,
Y is Y2 - Y1,
Z is Z2 - Z1.

normalize(Block, Norm) :-
minimum(Block, Min),
maplist(subtract(Min), Block, Translated),
sort(Translated, Norm).

all_rotations(Block, Rotations) :-
findall(Norm, (rotation(Block, Rot), normalize(Rot, Norm)), All),
sort(All, Rotations).

N > 0, N1 is N - 1,

close(Fd).

% Solver

solve(Model, Solution) :-
findall(B, block(B), Blocks),
solve(Blocks, Model, Solution).

solve([], [], []).
solve([Block|Rest], Model, [Transformed|Solution]) :-
all_rotations(Block, Rotations),
member(Rot, Rotations),
translate(Rot, Model, Transformed),
remove(Transformed, Model, Model1),
solve(Rest, Model1, Solution).

translate([Pos|Rest], Model, [Pos1|Trans]) :-
member(Pos1, Model),
subtract(Pos, Pos1, Delta),

X is X2 + X1,
Y is Y2 + Y1,
Z is Z2 + Z1.

remove([], Y, Y).
remove([X|Xs], Y, Z) :- select(X, Y, Y1), remove(Xs, Y1, Z).

% Output

writeVTK(Solution, Filename) :-
open(Filename, write, Fd),
write(Fd, '# vtk DataFile Version 2.0'), nl(Fd),
write(Fd, 'Block by Block solution'), nl(Fd),
write(Fd, 'ASCII'), nl(Fd),
write(Fd, 'DATASET POLYDATA'), nl(Fd),
write(Fd, 'POINTS 27 int'), nl(Fd),
write_points(Fd, Solution),
write(Fd, 'POINT_DATA 27'), nl(Fd),
write(Fd, 'SCALARS group int 1'), nl(Fd),
write(Fd, 'LOOKUP_TABLE default'), nl(Fd),
write(Fd, '1 1 1 1 2 2 2 2 3 3 3 3 4 4 4 4 5 5 5 5 6 6 6 6 7 7 7'), nl(Fd),
close(Fd).

write_points(_, []).
write_points(Fd, [[]|Rest]) :- write_points(Fd, Rest).
write_points(Fd, [[p(X,Y,Z)|Points]|Rest]) :-
write(Fd, X), write(Fd, ' '),
write(Fd, Y), write(Fd, ' '),
write(Fd, Z), nl(Fd),
write_points(Fd, [Points|Rest]).
```