I have a first solution Andre Dauphiné

enter code hereClearAll["Global`*"]
Needs["NDSolve`FEM`"]
carto = DiscretizeGraphics[
   CountryData["France", {"Polygon", "Mercator"}]];
bmesh = ToBoundaryMesh[carto, "MaxBoundaryCellMeasure" -> 25, 
   AccuracyGoal -> 1];
mesh = ToElementMesh[bmesh, MaxCellMeasure -> 5, 
  "MaxBoundaryCellMeasure" -> 25]
mesh["Wireframe"]
op = -Laplacian[u[x, y], {x, y}] - 20;
usol = NDSolveValue[{op == 1, DirichletCondition[u[x, y] == 0, True]},
   u, {x, y} \[Element] mesh]
Plot3D[usol[x, y], {x, y} \[Element] mesh]
Show[ContourPlot[usol[x, y], {x, y} \[Element] mesh, 
  ColorFunction -> "Temperature"], bmesh["Wireframe"]]