In a convolutional neural network (CNN), when convolving the image, is the operation used the dot product or the sum of element-wise multiplication? The example below is taken from the lectures in deeplearning.ai shows that the result is the sum of the element-by-element product (or "element-wise multiplication".  The red numbers represent the weights in the filter: 
$(1*1)+(1*0)+(1*1)+(0*0)+(1*1)+(1*0)+(0*1)+(0*0)+(1*1) 
= 1+0+1+0+1+0+0+0+1 
= 4 $

HOWEVER, most resources say that it's the dot product that's used:   

"…we can re-express the output of the neuron as , where is the bias
  term. In other words, we can compute the output by y=f(x*w) where b is
  the bias term. In other words, we can compute the output by 
  performing the dot product of the input and weight vectors, adding in
  the bias term to produce the logit, and then applying the
  transformation function."

Buduma, Nikhil; Locascio, Nicholas. Fundamentals of Deep Learning: Designing Next-Generation Machine Intelligence Algorithms (p. 8). O'Reilly Media. Kindle Edition. 

"We take the 5*5*3 filter and slide it over the complete image and
  along the way take the dot product between the filter and chunks of
  the input image.  For every dot product taken, the result is a
  scalar."

https://medium.com/technologymadeeasy/the-best-explanation-of-convolutional-neural-networks-on-the-internet-fbb8b1ad5df8

"Each neuron receives some inputs, performs a dot product and
  optionally follows it with a non-linearity."

http://cs231n.github.io/convolutional-networks/

"The result of a convolution is now equivalent to performing one large
  matrix multiply np.dot(W_row, X_col), which evaluates the dot product
  between every filter and every receptive field location."

http://cs231n.github.io/convolutional-networks/
However, when I research how to compute the dot product of matrics, it seems that the dot product is not the same as summing the element-by-element multiplication. What operation is actually used (element-by-element multiplication or the dot product?) and what is the primary difference? 
 A: I believe the key is that when the filter is convolving some part of the image (the "receptive field") each number in the filter (i.e. each weight) is first flattened into vector format.  Likewise, the pixels of the image are also flattened into vector format.  THEN, the dot product is calculated.  Which is the exact same thing as finding the sum of the element-by-element (element-wise) multiplication. 
Of course, these  flattened vectors can also be combined in a matrix format, as the below image shows.  In this case then true matrix multiplication can be used, but it's important to note that the flattening the image pixels from each convolution and also the weights filter is the precursor.  

image credit: TensorFlow and Deep Learning without a PhD, Part 1 (Google Cloud Next '17)
A: The operation is called convolution which involves a sum of element by element multiplication, which in turn is the same as a dot product on multidimensional matrices which ML folks call tensors. If you write it as a loop, it'll look like this pseudo Python code:
for i in 0:3:
  for j in 0:3:
    s = 0
    for k in 0:3:
      for l in 0:3:
        s += A[i+k,j+l]*C[k,l]
    Z[i,j] = s

Here A is your 5x5 input matrix, C is 3x3 filter and Z is 3x3 output matrix.
The subtle difference with a dot product is that usually a dot product is on the entire vectors, while in convolution you do dot product on the moving subset (window) of the input matrix, you could write it as follows to replace the innermost two nested loops in the code above:
Z[i,j] = dot(A[i:i+2,j:j+2],C)

A: Any given layer in a CNN has typically 3 dimensions (we'll call them height, width, depth). The convolution will produce a new layer with a new (or same) height, width and depth. The operation however is performed differently on the height/width and differently on the depth and this is what I think causes confusion.
Let's first see how the convolution operation on the height and width of the input matrix. This case is performed exactly as depicted in your image and is most certainly an element-wise multiplication of the two matrices.
In theory:
Two-dimensional (discrete) convolutions are calculated by the formula below:
$$C \left[ m, n \right] = \sum_u \sum_υ A \left[ m + u, n + υ\right] \cdot B \left[ u, υ \right]$$
As you can see each element of $C$ is calculated as the sum of the products of a single element of $A$ with a single element of $B$. This means that each element of $C$ is computed from the sum of the element-wise multiplication of $A$ and $B$.
In practice:
You could test the above example with any number of packages (I'll use scipy):
import numpy as np
from scipy.signal import convolve2d

A = np.array([[1,1,1,0,0],[0,1,1,1,0],[0,0,1,1,1],[0,0,1,1,0],[0,1,1,0,0]])
B = np.array([[1,0,1],[0,1,0],[1,0,1]])
C = convolve2d(A, B, 'valid')
print(C)

The code above will produce:
[[4 3 4]
 [2 4 3]
 [2 3 4]]

Now, the convolution operation on the depth of the input can actually be considered as a dot product as each element of the same height/width is multiplied with the same weight and they are summed together. This is most evident in the case of 1x1 convolutions (typically used to manipulate the depth of a layer without changing it's dimensions). This, however, is not part of a 2D convolution (from a mathematical viewpoint) but something convolutional layers do in CNNs.
Notes:
1: That being said I think most of the sources you provided have misleading explanations to say the least and are not correct. I wasn't aware so many sources have this operation (which is the most essential operation in CNNs) wrong. I guess it has something to do with the fact that convolutions sum the product between scalars and the product between two scalars is also called a dot product.
2: I think that the first reference refers to a Fully Connected layer  instead of a Convolutional layer. If that is the case, a FC layer does perform the dot product as stated. I don't have the rest of the context to confirm this.
tl;dr The image you provided is 100% correct on how the operation is performed, however this is not the full picture. CNN layers have 3 dimensions, two of which are handled as depicted. My suggestion would be to check up on how convolutional layers handle the depth of the input (the simplest case you could see are 1x1 convolutions).
A: Both element-wise and dot product interpretations are correct. When you convolve two tensors, X of shape (h, w, d) and Y of shape (h, w, d), you're doing element-wise multiplication. It's, however, the same as the dot product of X and Y transpose. You can expand the math equation, the shapes and subscripts match.
A: "What operation is actually used (element-by-element multiplication or the dot product?) and what is the primary difference?"
The actual operation is called a convolution defined as:
$$ J(x,y) = K * I = \sum_{n,m}K(n,m)I(x-n,y-m)  $$
Where $J$ is the convolved signal, $K$ is the kernel, $I$ is the input signal, and $n,m$ are the kernel indexes. We see by the input signal indexing $I(x-n,y-m)$ that both input signal axes are "flipped". As an aside, the kernel (a.k.a. filter) could be flipped instead.
Back to your question:
From the Numpy docs: the dot product numpy.dot "Returns the dot product of a and b. If a and b are both scalars or both 1-D arrays then a scalar is returned; otherwise an array is returned. (...)"
While the sum of the element-wise multiplication returns a scalar.
Making use of a previous reply - if we look at the two regions in the example:
import numpy as np
from scipy.signal import convolve2d

A = np.array([[1,1,1],
              [0,1,1],
              [0,0,1]])

B = np.array([[1,0,1],
              [0,1,0],
              [1,0,1]])

Dot product (returns ndarray):
print(np.dot(A,B))

[[2 1 2]
 [1 1 1]
 [1 0 1]]

Sum of element-wise multiplication (returns scalar);
print(np.sum(np.multiply(A,B)))

4

And that is the primary difference.
For good measure, a convolution:
D = convolve2d(A, B, 'valid')
print(D)

[[4]]

which in this case is equal to sum of element-wise multiplication of image patch and filter.
A convolution can be equivalent to sum of element-wise multiplication if the filter is symmetric, which is the case in deeplearning.ai example given. This will not be the case for asymmetric filters:
B2 = np.array([[1,0,0],
               [0,1,0],
               [1,0,1]])

print(np.sum(np.multiply(A,B2)))
print(convolve2d(A, B2, 'valid')) 

3
[[4]]

where sum of element-wise multiplication returns 3 and convolution returns 4, albeit as numpy.ndarray data type.
Note the discussion on flattened arrays - TensorFlow and Deep Learning without a PhD, Part 1 (Google Cloud Next '17) - applies to fully-connected, not convolutional layers.
2020 are you still reading? :D
