Index notation to matrix notation

Question

Assuming I have this equation in index notation and I am trying to convert it to a matrix notation: $$ A_{ij} = B_i {(C_j + D_{ij})}^\top $$ where $B, C \in \mathbb{R}^{N \times X}$ and $D \in \mathbb{R}^{N \times N \times X}$. Then I believe $A \in \mathbb{R}^{N \times N}$

If it wasn't for $D$, the answer would be pretty straightforward but the summation is really throwing me off here.

Here is an example of what I'm doing with Python code:

import numpy as np
A = np.zeros((2, 2))
B = np.random.randint(0, 4, (2, 3))
C = np.random.randint(0, 4, (2, 3))
D = np.random.randint(0, 4, (2, 2, 3))
for i in range(B.shape[0]):
  for j in range(C.shape[0]):
    A[i][j] = np.dot(B[i], C[j] + D[i][j])
B, C, D, A
(array([[2, 1, 3],
        [0, 0, 2]]),
 array([[2, 3, 0],
        [2, 0, 0]]),
 array([[[3, 3, 0],
         [0, 3, 1]],
 
        [[3, 1, 2],
         [0, 3, 0]]]),
 array([[16., 10.],
        [ 4.,  0.]]))

Is there a way to convert the above equation to a matrix notation?

Thanks!

EDIT: So I was able to come up with this code using numpy vectorization and it seems to work all right. However, I am not sure if it is mathematically sound or how to even translate it into proper mathematical notation.

np.diagonal(np.dot(B, D.transpose(0, 2, 1)), axis1=0, axis2=1).T + np.dot(B, C.T)