On a related question, namely how to simulate an iid random sample $(X_1,\ldots,X_n)\sim f$, under the constraint that its sum $X_1+\cdots+X_n$ is fixed to an arbitrary value $s_0$, I showed that the density of that sample is given by
$$(X_1,\ldots,X_n)\sim f(x_1)\cdots f(x_{n-1})f(s_0-x_1-\cdots-x_n)\mathbb{I}_{s_0-x_1-\cdots-x_{n-1}}(x_n)$$
Note: Under the constraint the sample is no longer independent but the
observations are identically distributed, despite the apparent
asymmetry in the above.
For instance, a Uniform sample with fixed average $0.5$ would have the joint density
$$(U_1,\ldots,U_N)\sim \prod_{i=1}^{N-1}\mathbb{I}_{(0,1)}(u_i)\mathbb{I}_{1/2}(u_1+\cdots+u_N)$$
which can simulated by Gibbs sampling or another MCMC algorithm like RWMH. When using Gibbs sampling, the $N-1$ first coordinates of the sample can be simulated one at a time with
$$U_i|U_{-i}\sim \mathbb{I}_{(0,1)}(u_i) \mathbb{I}_{(0,1)}(s_0-u_1-\cdots-u_i-\cdots-u_{n-1})=\mathbb{I}_{(\max\{0,s_0-1-\sum_{j\ne i,n}u_j\},\min\{1,s_0-\sum_{j\ne i,n}u_j\})}(u_i)$$
An R implementation looks like this code:
n=3;T=1e4
s0=.5 #fixed average
sampl=matrix(s0,T,n)
for (t in 2:T){
sampl[t,]=sampl[t-1,]
for (i in 1:(n-1)){
sampl[t,i]=runif(1,
min=max(0,n*s0-sum(sampl[t,c(-i,-n)])-1),
max=min(1,n*s0-sum(sampl[t,c(-i,-n)])))
sampl[t,n]=n*s0-sum(sampl[t,-n])}}
with the following marginals on the $U_i$'s:

Note: One can easily modify the above R code to impose an average of $s_0=0.05$ or $s_0=0.975$ on the sample.