Building on the trick of kguler and on the @user21 recommendations I ends up with the following approach. 

For many reason my actual interest is on first-order meshes and on a `MeshRegion`-type output. However, I used the some services of FEM toolbox and particularly of `NumericalRegion`.

 - Its `"BoundaryFunction"` property with `FindRoot` is useful to "improve" `RegionPlot` output. I don't know if we can easily build this compiled function directly from a generic `Region`.

 - For some misterious (for me) reason, the output of `RegionPlot[numericalRegion["Predicates"], ...]` appears more accurate than `RegionPlot[symbolicRegion, ...]`. This turns out to be very important to get a valid `MeshRegion` with the actual implementation because...

In the `RegionPlot` output improvement phase boundary vertices are moved, and a boundary vertex can cross a mesh line. So, at some times we need to reorder the vertices on a boundary, otherwise we can get an invalid / self-crossing boundary. This becomes a problem for fine meshes (in the following samples for `Mesh` option >= 50. A more accurate `RegionPlot` output is helpful to reduce the need of reordering.

I'd appreciate any suggestion to detect and fix this problem and to improve the whole implementation.


    << NDSolve`FEM`
    
    Options[ConstrainedRegionPlotMeshGenerator] = {Mesh -> Automatic, PlotPoints -> Automatic, MaxRecursion -> Automatic};
    ConstrainedRegionPlotMeshGenerator[region_NumericalRegion, opts : OptionsPattern[]] := 
     ConstrainedRegionPlotMeshGeneratorCore[region, opts] // ToBoundaryMesh
    ConstrainedRegionPlotMeshGenerator[region_ /; ConstantRegionQ[region],
       opts : OptionsPattern[]] := 
     ConstrainedRegionPlotMeshGeneratorCore[ToNumericalRegion@region, opts]
    
    ConstrainedRegionPlotMeshGeneratorCore[region_NumericalRegion, 
       opts : OptionsPattern[]] :=
      Module[{g, mptsx, mptsy, mpts, blns, mptq, grids, xg, yg, \[Delta], 
        pts, ptsx, ptsy, x, y},
       
       (* draw region with RegionPlot and extract the GraphicsComplex *)
       g = First@RegionPlot[
          (*region["SymbolicRegion"],*)
          region["Predicates"],
          Evaluate[
           Sequence @@ 
            MapThread[
             Prepend, {region["Bounds"], region["PredicateVariables"]}]],
          Frame -> False, PlotStyle -> None, BoundaryStyle -> Red, 
          MeshStyle -> Blue,
          Evaluate@FilterRules[{opts}, {Mesh, PlotPoints, MaxRecursion}]];
       
       (* get mesh lines endpoints and boundary lines *) 
       {mptsx, mptsy} = 
        Union @@@ 
         Cases[g, {Blue, lines___Line} :> {lines}, \[Infinity]][[All, All,
            1, {1, -1}]];
       If[Intersection[mptsx, mptsy] =!= {}, Return[$Failed]]; (* 
       Unsupported at present *)
       mpts = DeleteDuplicates@Join[mptsx, mptsy];
       blns = 
        Cases[g, {Directive[___, Red, ___], lines___Line} :> 
          lines, \[Infinity]];
       
       (* build closed boundary lines selecting boundary lines vertices \
    that are also on mesh lines *)
       mptq = AssociationThread[mpts, Range@Length@mpts];
       blns = Map[DeleteMissing[mptq /@ #] &, blns, {2}];
       blns = 
        Map[If[Last@# == First@#, #, Append[#, First@#]] &, blns, {2}];
       
       (* improve  boundary vertices position *)
       \[Delta] = region["BoundaryFunction"];
       ptsx = {#[[1]], 
           y /. FindRoot[\[Delta][{#[[1]], y}], {y, #[[2]]}]} & /@ 
         g[[1, mptsx]];
       ptsy = {x /. 
            FindRoot[\[Delta][{x, #[[2]]}], {x, #[[1]]}], #[[2]]} & /@ 
         g[[1, mptsy]];
       pts = Join[ptsx, ptsy];
       
       (* build result*)
       BoundaryMeshRegion[pts, Sequence @@ blns]
       ];



For example, with a uniform mesh, and a `MeshRegion` output:

    ConstrainedRegionPlotMeshGenerator[\[CapitalOmega], Mesh -> {5, 10}]

![Mathematica graphics](https://i.sstatic.net/jzJM5.png)


With an explicit, non-uniform grid, and a `MeshRegion` output:

    Show[
     ConstrainedRegionPlotMeshGenerator[\[CapitalOmega], Mesh -> grids],
     GridLines -> grids, GridLinesStyle -> LightGray, 
     Method -> {"GridLinesInFront" -> True}
     ]

![Mathematica graphics](https://i.sstatic.net/8uxCs.png)


With an explicit, non-uniform grid, and an `ElementMesh` output:

    mesh = ToBoundaryMesh[\[CapitalOmega], 
      "BoundaryMeshGenerator" -> {ConstrainedRegionPlotMeshGenerator, 
        Mesh -> grids}]
    
    Show[
     mesh["Wireframe"],
     mesh["Wireframe"["MeshElement" -> "PointElements"]],
     GridLines -> grids, GridLinesStyle -> LightGray
     ]

![Mathematica graphics](https://i.sstatic.net/whoEI.png)