HPLC: peak area vs concentration
However, in my opinion, there should be a proportionality between the concentration and the peak intensity, not the peak area.


There is a proportionality between both peak area vs. concentration and peak height vs. concentration. 


Peak height is proportional to the instantaneous amount of analyte that is transiting the detector.  
Peak area is proportional to the sum of all of analyte moleucles that have transited the detector.  
From 1 and 2 you might be able to infer the relationship betweeen peak height $h$ and area $A$: $$A = \int{h(t)\;dt}$$
People are usually interested in the total amount of substance injected into the column.  (If they know the injection volume, they can calculate the concentration from this value.)  The total amount of substance would be calculated from the peak area.
The maximal peak height is proportional to the peak area, but only if the peak "shape" is constant.  Here are some scenarios where peak shape will not be constant:


You want to compare injection #2 you made on your column two years ago to injection #2000 that you made recently.  Due to column degradation, the recent injections have much more noticeable peak tailing.  Because the longer tail leads to wider, asymmetric peaks, there will be less total height at the peak maximum relative to "better", symmetric peaks.
You change the flow rate of the HPLC.  All peaks are narrower, and thus, higher.
The scan rate of your detector changes, and the detector reports "counts" or "intensity" rather than counts per time.  This is actually very common for mass spectrometers, but less common for absorbance detectors.  If you scan twice as fast in MS, your peak heights go down by approximately twofold, but you have data points twice as often, so the peak area is relatively unchanged.  
Your detector undersamples the peak.  Say your chromatography is nice and the "real" peak shape is nicely gaussian, but is ~ten seconds wide and your detector only gives you a data point every two seconds.  Say that due to very small random drifts in retention time, on some injections one of the data points coincides exactly with the maximum of the "real" peak, but in others, it is a little bit off.  The peak heights in this scenario will vary considerably more than the peak areas. (If the recorded maximum is off-center from the true maximum, there will be two data points that are higher, i.e. "closer" to the maximum than if the data maximum coincides with the true maximum, so integration will partially correct this error.)  Essentially, the peak maximum is a single-point sample from the gaussian distribution of the peak, while the area is a several-point sample from the gaussian distribution, so it has better sampling properties.




  Is there any mathematical relation between the peak area and the concentration?


Yes, there is, but it depends on the peak shape.  For perfectly gaussian peaks, $$h_{max} = \frac{A}{\sigma \sqrt{\tau}}$$

where $\tau$ is $2 \pi$ and $\sigma$ is the width of the peak, which is related to the full-width at half-maximum of the peak by $\mathrm{FWHM} = 2 \sigma \sqrt{2 \ln 2}$.

However, for non-gaussian peaks, this relationship does not hold.