General Criteria for the Existence or Non-existence of Solutions to a PDE What are the conditions[General Criteria] for the existence or non existence of the solutions to a PDE[Elliptic type] subject to given boundary conditions?
A specific Example:
Let's consider the reduced wave equation:
$$\frac{\partial^2 u}{\partial x^2}+\frac{\partial^2 u}{\partial y^2}=-k^2u$$
In the simplest case if u vanishes on the square:$x=0,x=a,y=0,y=a$ the solution is:
$$u=\sin\frac{m\pi x}{a}\sin\frac{n\pi\ y}{a}$$
where,
$k^2=(m^2+n^2)\pi^2/a^2$.
Can we use this result to solve:
$$\frac{\partial^2 u}{\partial x^2}+\frac{\partial^2 u}{\partial y^2}=f(u)$$
One may to expand f(u) in terms of the eigenfunctions of the operator:
$$\frac{\partial^2 }{\partial x^2}+\frac{\partial^2 }{\partial y^2}$$
How should one go about the job if the boundary conditions on the same square are changed to a different continuous and differentiable functions of the type:
Example:
For x=0 $u=f_1(x,y)$
For x=a $u=f_2(x,y)$
For y=0 $u=f_3(x,y)$
For x=a $u=f_4(x,y)$
 A: The answer to the general question you asked

What are the conditions[General Criteria] for the existence or non existence of the solutions to a PDE[Elliptic type] subject to given boundary conditions?

Is "there are none". To be a bit more verbose: the answer is that "it depends on what form of equation" you are allowed to consider and "what you admit as solutions". 
Let me give one extremely general example. Consider the family of partial differential equations $F(\partial^2 u) = 0$, where $F$ is some smooth function defined on the space of $N\times N$ matrices. The equation is said to be elliptic if $\frac{\partial}{\partial X_{ij}} F(X_{ij})$ is elliptic (as a matrix, evaluated at all $X$). For these equations, 


*

*In dimension $N = 2$, Nirenberg showed that the Dirichlet problem admits classical solutions. 

*In arbitrary dimensions, viscosity solutions exist (using the Perron method). 

*If $F$ is convex, viscosity solutions are in fact classical (Krylov-Evans regularity).

*If $F$ is not convex, there can be cases where a viscosity solution exists but is not classical, which by the uniqueness properties of viscosity solutions, imply that such equations do not admit classical solutions. 


If you restrict to only quasilinear instead of fully nonlinear equations, there are a wealth of existence results for the Dirichlet problem. In fact, in section 15.5 of their book, Gilbarg and Trudinger wrote that

It is not feasible to present here a comprehensive account of existence theorems for the classical Dirichlet problem that follow by combination of the results of Chapters 10, 13, 14, and 15. 

... and they proceed to list 10 different existence theorems applicable in different contexts. What's important to note, however, is that there is a fairly strong non-existence theorem that applies even in the case of the semilinear problem. I will state here a special case
Theorem (c.f. Theorem 14.11 in Gilbarg-Trudinger) Let $\Omega\subset\mathbb{R}^n$ be a bounded domain such that the largest ball contained in $\Omega$ has radius $R$, and consider the Dirichlet problem for 
$$ \triangle u(x) + b(x,u, \nabla u) = 0~. $$
If we assume that there is a positive constant $\theta$ such that for every $p\in \mathbb{R}^n$, $z \in \mathbb{R}$, and $x\in\mathbb{R}^n$ you have that 
$$ |b(x,z,p)| \geq \frac{n |p|}{R} + |p|^{2+\theta} $$
then one can find a smooth boundary data $\phi$ for which the Dirichlet problem does not admit a solution. 
Note, however, that this theorem does not apply to your special case, since your $f(u)$ does not depend on the first derivative. 

Now, to address your specific problem: you outlined an approach that requires expanding $u$ in terms of the eigenfunctions of the Laplacian. This is roughly equivalent to taking the Fourier transform of the problem, and so you see that to directly apply the solutions to the reduced wave equation you need that $f(u)$ to be linear ($f(au+bv) = af(u) + bf(v)$), which puts strong restrictions on $f$. 
Now, for a semilinear problem, one can always convert a problem with nonzero boundary condition to one with zero boundary conditions: say you want to solve
$$ \triangle u = f(x,u) $$
with boundary condition $\phi$. First solve the homogeneous problem $\triangle w = 0$ with boundary condition $\phi$. Then write the ansatz that $u = v + w$ where $v$ vanishes on the boundary. Then $v$ solves
$$ \triangle v = \tilde{f}(x,v) := f(x,v+w) $$
In general, however, for a problem like the one you proposed $\triangle u = f(u)$, you are often better off demonstrating existence of a solution using a variational method (if possible) or a fixed-point/iteration method. (If you are really lucky you can use Perron method, but for general nonlinearities $f(u)$ the comparison results for super and sub solutions may not be available.) 
A: In Response to the Specific problem in the Question:
On the evaluation of particular Integrals:
$$\frac{\partial^2 u}{\partial x^2}+\frac{\partial^2 u}{\partial y^2}=f(u)$$
We look for a particular soln: dependent on x only:$u=u(x)$
It represents a 2D surface whose gradient changes only along the x-direction but not along the y-direction.
We solve the ODE:
$$\frac{d^2 u}{dx^2}=f(u)$$
$$<=>pdp=f(u)du$$
Where $p=\frac{du}{dx}$
In several cases we may obtain an analytical soln for:
$$\int f(u)du$$
This particular integral may contain arbitrary constants which we may use to our best advantage in catering to some given boundary condition.
Similarly,we may start off with a solution of the form
$$u=u(y)$$
To obtain another class of particular solutions.
The functions u=u(x) provide 2D surfaces with slope changing in the x-direction while v=v(y) provide surfaces with slope changing in the y-direction. These surfaces will be criss-crossing each other. We may consider infinitesimally small surfaces of the two types and sew them into a smooth, continuous curved surface which has gradient changing both in the x and the y directions at different points
PDE:$$\frac{\partial^2 u}{\partial x^2}+\frac{\partial^2 u}{\partial y^2}=\cot(u)$$
 Trial Soln: $u=2x^5+3y^7+\lambda x^2$ -------------(1).
Using this relation in the PDE we have $40x^3 +126y^5 +2λ=cot(u)$ -------(2).
But (1) and (2) represent different surfaces. The PDE is valid on the first surface at the intersection of the two surfaces for an arbitrary λ . But we have our handle on λ . We may change it by infinitesimal amounts to generate a solution surface.We may use arbitrary trial solutions of the type:$$u=Af_1(x)+Bf_2(y)+\lambda x^2$$.
The constants A and B are to be used for fitting in boundary conditions while $\lambda$ is to be used for generating infinitesimal surface strips.
We again consider the PDE with f(u)=cotu.
PDE:$\frac{\partial^2 u}{\partial x^2}+\frac{\partial^2 u}{\partial y^2}=cot(u)$. 
Soln: $\int_a ^u \frac {dt}{\sqrt{ln\mid K\sin(t)\mid}}=Ax+By+C$----(1),
where $A^2+B^2=2$.
By differentiating (1) we have,
$\frac{\partial u}{\partial x}=A\sqrt{ln\mid K\sin(u) \mid}$
$\frac{\partial u}{\partial y}=B\sqrt{ln\mid K\sin(u) \mid}$
$\frac{\partial^2 u}{\partial x^2}=(A^2/2) \cot(u)$
$\frac{\partial^2 u}{\partial y^2}=(B^2/2) \cot(u)$
The solution indicated by relation (1) indeed satisfies the PDE.
A: Let’s solve the PDE:$$\frac{\partial^2 u}{\partial x^2}+ \frac{\partial^2 u}{\partial y^2}=\cot (u)$$  ------------(1) 
Subject to the constraint :$$(\frac{\partial u}{\partial x})^2+(\frac{\partial u}{\partial y})^2=1$$ ----------------- (2)
We write PDE (1) as: $$\frac{\partial^2 [(u+F(u))-F(u)]}{\partial x^2}+ \frac{\partial^2 [(u+F(u))-F(u)]}{\partial y^2}=\cot (u)$$ ----------- (3)
PDE (3) may be represented as the following two equations:
$$\frac{\partial^2 (u+F(u))}{\partial x^2}+ \frac{\partial^2 (u+F(u))}{\partial y^2}=0$$--------- (4) 
And
$$\frac{\partial^2 F(u)}{\partial x^2}+ \frac{\partial^2 F(u)}{\partial y^2}=-\cot (u)$$--------------- (5)
Now,
$$\frac{\partial F(u)}{\partial x}=\frac{dF}{du} \frac{\partial u}{\partial x}$$
$$\frac{\partial^2 F(u)}{\partial x^2}=\frac {d^2 F}{du^2} (\frac {\partial u}{\partial x})^2+\frac{dF}{du}\frac{\partial^2 u}{\partial x^2}$$
And
$$\frac{\partial F(u)}{\partial y}=\frac{dF}{du} \frac{\partial u}{\partial y}$$
$$\frac{\partial^2 F(u)}{\partial y^2}=\frac {d^2 F}{du^2} (\frac {\partial u}{\partial y})^2+\frac{dF}{du}\frac{\partial^2 u}{\partial y^2}$$
PDE (5) may be written as:
$$\frac {d^2 F}{du^2} [(\frac {\partial u}{\partial x})^2+(\frac {\partial u}{\partial y})^2]+\frac{dF}{du}[\frac{\partial^2 u}{\partial x^2}+\frac{\partial^2 u}{\partial y^2}]=-\cot(u)$$
Or,
$$\frac {d^2 F}{du^2}\times {1}+\frac{dF}{du}\times \cot(u) =-\cot(u)$$
Or,
$$\frac {d^2 F}{du^2}+\cot(u)\frac{dF}{du} =-\cot(u)$$
The above ODE gives us the form for F to be used in PDE (4) which is Laplace’s equation in two dimensions
We obtain: $$ u+F(u)=Soln \;of \;Laplace’s\; eqn.\; in \; two \;dimensions$$
[Differentiating the above relation twice wrt to x and y and adding the results we have,
$$F''(u)[(\frac{\partial u}{\partial x})^2+(\frac{\partial u}{\partial y})^2]+F'(u)[\frac{\partial^2 u}{\partial x^2}+\frac{\partial^2 u}{\partial y^2}]=-[\frac{\partial^2 u}{\partial x^2}+\frac{\partial^2 u}{\partial y^2}]$$
Again F(u) satisfies relation (5):$$F''(u)+\cot(u)F'(u)=-\cot(u)$$
 PDE's (1) and (2)  getting satisfied simultaneously is one of the options--the simplest one possibly]
In the final step the particular solutions for F and that for Laplace's equation are chosen in a manner such that equation (2) is satisfied.
