Generate three random numbers that sum to 1 in R I would like to generate three random numbers and then standardize them so that they add up to 1.
I would like to repeat this procedure so that in the long run the mode is .33 for each number. 
 A: If X1, X2, and X3 are i.i.d. Gamma(a) then {X1,X2,X3}/(X1+X2+X3) will be Dirichlet(a,a,a). 
If a>1 then the mode will be 1/3. The peak will be sharper for larger values of a.
A: The mode is a bit of a red herring. Here is a very simple solution to this problem that circumvent the need to define the mode precisely. I'm surprised it has not been proposed earlier. The constraint on the mode can be easily satisfied by drawing samples from a symmetric distribution and scaling them suitably:
$$(x_i,y_i,z_i)\sim\mathrm{i.i.d.}\;\mathcal{L}(\mu,\sigma)$$ 
$$(x_i^*,y_i^*,z_i^*)=\left(\frac{x_i}{x_i+y_i+z_i},\frac{y_i}{x_i+y_i+z_i},\frac{z_i}{x_i+y_i+z_i}\right)$$
where $\mathcal{L}(\mu,\sigma)$ is a symmetric distribution (so that the mean, the mode and the median are the same) and chosen such that the probability mass below 0 is 0.  For example, picking $\mathcal{L}(\mu,\sigma)$ to be $\mathrm{Beta}(2,2)$:
a1 <- matrix(rbeta(100*3,2,2), nc=3)
a1 <- sweep(a1, 1, rowSums(a1), FUN="/")
colMeans(a1)
# [1] 0.3342165 0.3341534 0.3316301

yielding the desired solution
sum(colMeans(a1))
# [1] 1

A: Here is an approximate numerical answer. It can easily be made more precise.
Let $\{U,V,W\} = {X,Y,Z}/(X+Y+Z)$, where $X,Y,Z$ are i.i.d. with a trapezoidal density on $[0,1]$:
$f(x)=1+a-2ax.$ $U,V,W$ will have identical marginals.
Given a numeric 'a', I used Mathematica to get the cdf of $U$:
F[u_] = Assuming[0 < u < 1, Simplify@Integrate[  
Boole[x < u(x+y+z)] f[x] f[y] f[z], {x,0,1},{y,0,1},{z,0,1}]

Differentiating $F$ twice, setting the result to zero, and solving the resulting 7th degree polynomial gave the mode. I used a binary search to refine the value of 'a'. I used exact arithmetic throughout, up to the point of solving the polynomial.
  a      mode

  1    .318182  
 7/8   .322065  
13/16  .327099  
25/32  .330465  
49/64  .332373  
97/128 .333376 <-- close enough?  
 3/4   .334221  
 1/2   .353738  
  0    .359187

A: It is still unclear whether the OP wants a solution with a mode of 0.33 or $1 \over 3$ or a mean with one of those two values. Without knowing the exact need, there are multiple possibilities. [1] and [4] below address the problem of getting a mean of ${1 \over 3},$ while [2] and [3] are for a mode of ${1 \over 3}$. 
[1] Generate $U_1, U_2, U_3$ as continuous uniform random variates on $[0,{2 \over 3}]$. Let $A = {{U_1} \over {U_1 + U_2 + U_3}},$ $B = {{U_2} \over {U_1 + U_2 + U_3}},$ and $C = {{U_3} \over {U_1 + U_2 + U_3}}.$
[2] Let $X=(1/3)*W_1 + (1/6),$ $Y = (1/3)*W_2 + (1/6),$ and $Z = 1 - X - Y,$ where $W_1$ and $W_2$ are continuous uniform on $[0,1].$ $Z$ has a different distribution than $X$ or $Y,$ but I think this will meet the original mode requirement. 
[3] In an effort to produce a more intuitive version, let $R_1$ be right triangular with left endpoint at zero and mode at ${1 \over 3} .$ Let $L_1$ be left triangular with mode at $1 \over 3$ and right endpoint at ${2 \over 3} .$  Then $R_1 + L_1$ has a unique mode at ${2 \over 3},$ and if we define $Q = 1 - (R_1 + L_1)$ then $Q$ has a unique mode at ${1 \over 3}.$ The pdf of $Q$ is given below in the comments. 
[4] Trying to be clever, simple, and elegant. Generate 2 independent uniform[0,1] realizations. These 2 points divide the interval from 0 to 1 into 3 pieces. Use the lengths of these 3 pieces as the desired variates. Note how this generalizes intuitively to any sum and any number of random variables. Each variate is identically distributed (another right triangular distribution). Every pairwise correlation is $-{1 \over 2}.$ However, like approach [1], the mean is ${1 \over 3},$ but the mode here is not ${1 \over 3}.$ As whuber noted, it is at zero. 
A: Analytically: Given a joint pdf for $X$,$Y$, and $Z$ $f_{X,Y,Z}(x,y,z)$, if they are iid, then $f_{X,Y,Z}(x,y,z)=f_X(x)f_Y(y)f_Z(z)$, where $f_X(x)=f_Y(y)=f_Z(z)$. You have to find the pdf $$f_{X,Y,Z}\left(\dfrac{x}{x+y+z}\right)$$ After differentiating and equal to zero, you'll find your mode. Obviosly, mode will depend on $f_{X,Y,Z}(x,y,z)$ and consequently on $f_X(x)$, $f_Y(y)$ and $f_Z(z)$.
Numerically: sample three random numbers from your preferred distribution, standardize them and save them in the $i^{th}$-row an Nx3 matrix. Repeat this procedure N times and plot the frequencies of each column.
The analytical solution is preferred instead of trying to demonstrate it from random samples in R, which would be just a numerical approximation.
A: Here is a simple approach. Generate $X \sim \mathrm{Unif} [0,{2 \over 3}].$ Let $$Y = \begin{cases} X+{1/3} \ , & \text{if} \ X \le {1/3} \\ X-{1/3} \ , & \text{if} \ {1/3} \lt X \le {2/3}  \end{cases}$$
Let $Z = 1 - X - Y.$ 
Then it's not too hard to show that $X,Y,$ and $Z$ are identically distributed with $\mathrm{Unif} [0,{2 \over 3}]$ distributions. So each has mean and median of ${1 \over 3}$ and has a mode there as well. 
Additionally, all pairwise correlations are $ -{1 \over   2},$ and only one call to a uniform random generator is needed to get the three variates.  
