I want to generate a small random perturbation on a flat surface defined on a square domain, e.g. 1 + 0.05 rand(x,y), which can be used as initial condition with periodic boundary conditions in NDSolve. In other word, rand(x,y) is a special pseudo-random function distributed in the interval (-1,1) which has periodicity at the sides of the square-domain [0,L]*[0,L], that is, rand(0,y)=rand(L,y) and rand(x,0)=rand(x,L). In fact, this question is very similar to this question. However, after I went through this post in details, all comment, voted answers and even "accepted" answer as well as the answer given by the author, I found this question is still unsolved. Based on acl's comment and Vitaliy's answer, I have the following working example with two different initial conditions.

Off[NDSolve::mxsst];
L = 4*\[Pi];
c = -(1/20);
tmax = 1000;
ini1 = Interpolation@Flatten[Table[{{x, y}, 1 + c*RandomReal[{-1, 1}]}, {x, 0, L + 1}, {y, 0, L + 1}], 1]; (*initial condition given by acl*)
cond3D = First[h /. NDSolve[{D[h[x, y, t], t] + 
Div[h[x, y, t]^3*Grad[Laplacian[h[x, y, t], {x, y}], {x, y}], {x, y}] + 
Div[h[x, y, t]^3*Grad[h[x, y, t], {x, y}], {x, y}] == 0,
h[0, y, t] == h[L, y, t],
h[x, 0, t] == h[x, L, t],
h[x, y, 0] == ini1[x, y]
 },
h, {x, 0, L}, {y, 0, L}, {t, 0, tmax},
Method -> {"BDF", "MaxDifferenceOrder" -> 1}, 
MaxStepFraction -> 1/50]]

I have the following error:

NDSolve::ibcinc: Warning: boundary and initial conditions are inconsistent.

NOTE: When you try the second initial condition please save all the important files, because it is very easy to make your machine run away! I am sooo... sorry!

If I use the other initial condition given by @Vitaliy just with a slight modification:

ini2 = 1 + c*BSplineFunction[RandomReal[{-1, 1}, {30, 30, 1}], SplineClosed -> True];

where, because I want a distributed in the interval (-1,1) so {-1,1}; and because I want a small perturbation so a small parameter c. Just replace the initial condition as

h[x, y, 0] == ini2[x, y]

The NDSlove then eats up the memory and spit a lot a error related to data.

I realized that Vitaliy's answer is wonderful in SplineClosed -> True that makes I can use it with periodic boundary conditions in NDSolve. But another confused issue is when I try, ini2[9, 0] === ini2[9, L] // FullSimplify using my ini2, MMA gives False.I know @Vitaliy produced a disturbed surface in the domain[0, 1]* [0, 1], but how to get such a surface in a square domain, [0, L]*[0, L].

How can I generate a random disturbance consisting with periodic boundary condition meanwhile. If anyone can give such an initial condition that has something like the SplineClosed->True feature that can be fed in NDSolve, that would be truly awesome! I admit I am being very silly when I try to use RandomReal[] in this way.