Regarding your plots, I assume that fitting `aov()` with the `Error()` function won't work because you will get more than one error stratum.
You can see this with `print(ex.aov)`. Now you could use the `proj()` function which will give you the residuals
for each error stratum in a way that allows you to extract them more easily.
I found some information [here](http://stackoverflow.com/questions/26169153/how-to-get-residuals-from-repeated-measures-anova-model-in-r).

    print(ex.aov)

    ex.aov.proj <- proj(ex.aov)

    # Check for normality by using the 5th error stratum as an example since there is no `Error:Within` stratum
    qqnorm(ex.aov.proj[[5]][, "Residuals"])
    # Check for heteroscedasticity by using the 5th error stratum as an example since there is no `Error:Within` stratum
    plot(ex.aov.proj[[5]][, "Residuals"])

However, this will also lead into plots which I cannot fully interpret but I think it has to do with the fact that you don't
have multiple measurements for all factor combinations within `p` and your error term doesn't allow for an `Error:Within` stratum. But I am not sure. Does your real dataset have the exact
same structure (inlcuding number of observations)? Hopefully someone else can clarify.

**My alternative suggestion:**

First, I changed your dataset slightly and set a seed to make it reproducible (might be handy for some problems you have in the future):
 
    # Set seed to make it reproducible
    set.seed(12)

    # I changed the names of your variables to make them easier to remember
    # I also deleted a few nested `rep()` commands. Have a look at the `each=` argument.
    subj <- sort(factor(rep(1:20,8)))
    x1 <- rep(c('A','B'),80)
    x2 <- rep(c('A','B'),20,each=2)
    x3 <- rep(c('A','B'),10, each=4)
    outcome <- rnorm(80,10,2)

    d3 <- data.frame(outcome,subj,x1,x2,x3)

Second, I used a linear mixed-effects model instead since you have repeated measures and hence a random term you can use:

    require(lme4)
    # I specified `subj` as random term to account for the repeated measurements on subject.
    m.lmer<-lmer(outcome ~ x1*x2*x3 + (1|subj), data = d3)
    summary(m.lmer)

    # Check for heteroscedasticity
    plot(m.lmer)

[![enter image description here][1]][1]

    # or
    boxplot(residuals(m.lmer) ~ d3$x1 + d3$x2 + d3$x3)

[![enter image description here][2]][2]

    # Check for normality
    qqnorm(residuals(m.lmer))

[![enter image description here][3]][3]

Using the `afex` package you can also get the fixed effects in ANOVA table format (you can also use the `Anova()` function from the `car` package as another option):

    require(afex)
    mixed(outcome ~ x1*x2*x3 + (1|subj), data = d3, method="LRT")

    Fitting 8 (g)lmer() models:
    [........]
        Effect df    Chisq p.value
    1       x1  1     0.04     .84
    2       x2  1     2.53     .11
    3       x3  1  7.68 **    .006
    4    x1:x2  1  8.34 **    .004
    5    x1:x3  1 10.51 **    .001
    6    x2:x3  1     0.31     .58
    7 x1:x2:x3  1     0.12     .73

Check `?mixed` for the various options you can choose. Also  regarding mixed models, there is a lot of information here on Cross Validated.


  [1]: https://i.sstatic.net/i1X6E.png
  [2]: https://i.sstatic.net/k58cP.png
  [3]: https://i.sstatic.net/R6PRH.png