Constrained block placement puzzle

I'm trying to design a mathematical model for an algorithm I need for my program.

The problem

I have a 3d grid in which I have a shape $$S$$ and I need to fill this shape with block types $$B$$ in a way that every block type $$b \in B$$ satisfy its own 6 neighbor constraints. In other words, for every $$b$$, for each of the 6 directions $$D$$, I have a list of block types that are allowed to be placed nearby and I need to arrange these block types in order to fill the shape.

Initial steps

I've found a paper that describes a similar problem (The Jigsaw Puzzle) and propose a mathematical model for it: https://shaharkov.github.io/projects/GlobalPuzzles.pdf.

My first guess is: I need for every position of $$S$$ to know which $$b$$ to associate so the ideal thing would be to have a system of $$N$$ equations ($$N = \#S$$), with $$N$$ variables $$(b_{p_0}, ..., b_{p_{n-1}})$$ with $$b_{p_0}$$ being the block type to set for position 0.

However what resulted to be hard is writing the system itself:

$$$$\begin{cases} \sum_{d=1}^{6}r_d(b_{p_0}, b_{rel(p_0, d)})) = 6 & p_0 = (0, 0, 0)\\ .\\ .\\ .\\ \sum_{d=1}^{6}r_d(b_{p_{n-1}}, b_{rel(p_{n-1}, d)})) = 6 & p_{n-1} = (p_{{n-1}_x}, p_{{n-1}_y}, p_{{n-1}_z})\\ \end{cases}$$$$

I imagined to have a $$r_d$$ function (with $$d$$ being one of the 6 directions: front, back, left, right, top, bottom), that given two block types - one at position $$p$$ and the other at the neighbor position subject to $$d$$ - returns $$\{0, 1\}$$, $$0$$ if the two blocks doesn't match, $$1$$ if the two blocks match.

However it turned out it's really hard to think how to implement $$r_d$$: for a block type at a certain position $$b_p$$ I would need to mathematically encode that $$b_{rel(p, d)}$$ must be in the list of allowed blocks of $$b_p$$ for direction $$d$$ and therefore the function can return $$1$$, otherwise $$0$$.

Can anyone help me with this?