def assemble(self):
        ndim = self.mats[0].dimensions
        shape = self.mats[0].space.shape(True)
        self.solvers1D = []
        if ndim == 2:
            zi = np.ndindex((1, shape[1])) if self.naxes == 0 else np.ndindex((shape[0], 1))
            other_axis = (self.naxes+1) % 2
            for i in zi:
                sol = None
                for mat in self.mats:
                    sc = mat.scale[i] if mat.scale.shape[other_axis] > 1 else mat.scale[0, 0]
                    if sol:
                        sol += mat.mats[self.naxes]*sc
                    else:
                        sol = mat.mats[self.naxes]*sc
                self.solvers1D.append(Solver(sol))

        elif ndim == 3:
            s = [0, 0, 0]
            n0, n1 = np.setxor1d((0, 1, 2), self.naxes)
            for i in range(shape[n0]):
                self.solvers1D.append([])
                s[n0] = i
                for j in range(shape[n1]):
                    sol = None
                    s[n1] = j
                    for mat in self.mats:
                        sc = np.broadcast_to(mat.scale, shape)[tuple(s)]
                        if sol:
                            sol += mat.mats[self.naxes]*sc
                        else:
                            sol = mat.mats[self.naxes]*sc
                    self.solvers1D[-1].append(Solver(sol))