Define Implicit Region using Interpolating Function

Question

Suppose I have a set of points $\{(x_1,x_2) : g(x_1,x_2) > f(x_1,x_2)\}$, how would I turn this into a mesh that acts as a region for PDE solvers, or even to use as a plotting region?

I have tried using Implicit Region, but my function $f$ is an interpolating function and so is defined over a Disk region, which is not a square domain.

Hence if I run

omega = ImplicitRegion[g[x1,x2] > f[x1,x2], {{x1, -1, 1}, {x2, -1, 1}}]

and try to plot my region

Region[omega]

It doesn't give me anything, and trying to use this omega for use as a plotting domain / PDE domain doesn't work either. However, I can plot my region by running

set = f[x1,x2] > g[x1,x2]
RegionPlot[set, {x1, -1, 1}, {x2, -1, 1}]

I suppose this means RegionPlot is more robust in the sense that it doesn't care if my function $f$ is indeterminate at certain points.

I have also tried to change the values of $x_1,x_2$ in the implicit region definition to where my function $f$ is defined i.e. commands like

omega = ImplicitRegion[g[x1,x2] > f[x1,x2], {{x1,x2} \in Disk[0,0,r]}]

dont seem to work for ImplicitRegion.

Here is a minimal working example:

sol = NDSolve[{-Laplacian[u[x, y], {x, y}] == 1, 
    DirichletCondition[u[x, y] == 0, True]}, 
   u, {x, y} \[Element] Disk[]];
g[x1_, x2_] := x1^2 + x2^2;
f[x1_, x2_] := u[x1, x2] /. sol[[1]];
setfn[x1_, x2_] := g[x1, x2] > f[x1, x2];
omega = ImplicitRegion[setfn[x1, x2], {x1, x2}];
Region[omega]
Plot3D[{f[x1, x2], g[x1, x2]}, {x1, x2} \[Element] Disk[]]

My problem is that the Region[omega] part of the code gives me some weird Region[Embedding Dimension 2] output instead of actually plotting the region, whereas

RegionPlot[setfn[x1, x2], {x1, -1, 1}, {x2, -1, 1}]

plots the region, even though $f$ is undefined at points outside the disk.

I want to use such an omega as the domain to solve my PDE over, for example

 NDSolve[{-Laplacian[u[x, y], {x, y}] == 0, 
  DirichletCondition[u[x, y] == Sin[x y], True]}, u, {x, 
   y} \[Element] omega]

Also tried to extend my interpolating function with values that are irrelevant to my problem

f[x1_, x2_] := 
 Piecewise[{{u[x1, x2] /. sol[[1]], {x1, x2} \[Element] Disk[]}, {0, 
    True}}]

so it is defined everywhere. This doesn't seem to change it. Perhaps ImplicitRegion does not like interpolating functions at all?

Answer 1

For a subsequent FEM analysis you can use

Needs["NDSolve`FEM`"]
mesh = ToElementMesh[omega];
mesh["Wireframe"]

For other region related computations use:

DiscretizeRegion[omega]