Let's denote the enclosing box by $\varOmega$, its boundary by $\varGamma_1 = \partial \varOmega$, the polygonal line by $\varGamma_2$, and the probability distribution by $u$.
Let $a$ denote the diffusivity of the Brownian process
and let $D_i$ and $N_i$ denote the Dirichlet and Neumann operators (the latter with respect to $a$) of $\varGamma_i$, respectively, and let $\mathcal{H}^k$ denote the $k$-dimensional Hausdorff-measure.
For me it sounds as if $u$ had to satisfy the following linear PDE:
$$\begin{array}{rcll}
- \operatorname{div}(a \, \operatorname{grad} u)
&=
& - \frac{1}{\mathcal{H}^2(\varOmega)} \int_{\varGamma_2} (N_2 \, u) \, \operatorname{d} \mathcal{H}^1,
&\text{(Gauß' law: matter orginates from teleporting)}\\
N_1\,u &= &0,
&\text{(box is isolated)}
\\
D_2\,u &= &0,
&\text{(matter is immediately teleported away)}
\\
\int_\varOmega u \, \operatorname{d} \mathcal{H}^2 &= &1.
&\text{(gauging: $u$ has to be a probability density)}
\end{array}
$$
While the other three equations might be evident, here a short explanation why I think that the first equation has to hold for the steady state solution $u$ of the Brownian motion: The rate $\int_{\varGamma_1} N_2 \,u \, \operatorname{\mathcal{H}}^1$ of particles hitting $\varGamma_2$ has to be equal to the rate of particles appearing everywhere else. By Gauß' law, the latter is $\int_\varOmega ( - \operatorname{div}(A \, \operatorname{grad} u)) \, \operatorname{\mathcal{H}}^2$.
Actually, this should be solvable by the linear system for $(u,\lambda)$:
$$\begin{array}{rcll}
- \operatorname{div}(a \, \operatorname{grad} u) &= & \lambda,
\\
N_1\,u &= &0,
\\
D_2\,u &= &0,
\\
\int_\varOmega u \, \operatorname{d} \mathcal{H}^2 &= &1.
\end{array}
$$
The mass balance would automatically imply
$$\int_\varOmega \lambda \, \operatorname{d}\!x + \int_{\varGamma_2} (N_2 \, u) \, \operatorname{d}\!s = 0.$$
Here $\lambda \in \mathbb{R}$ works as the Lagrange multiplier for the probability conservation equation $\int_\varOmega u \, \operatorname{d} \mathcal{H}^2 = 1$.
Numerical solution
This PDE can be solved by means of finite elements. Due to the right-hand side of the first equation, this is maybe not possible with the high-level facilities in Mathematica (i.e., NDSolve) alone. But the low-level facilities should be able to provide us with the system matrix for this linear equation; solving it with LinearSolve is a standard task.
Let's start by defining the mesh on which to perform computations.
Needs["TriangleLink`"];
(* half of outer box's edge length*)
L = 5;
(* half of inner box's edge length*)
a = 1;
(* subdivision count for outer and inner box, respectlively*)
{m, n} = {100, 50};
h1 = N[L/m];
h2 = N[a/n];
Γ1 = DiscretizeRegion[
RegionBoundary@Rectangle[{-L, -L}, {L, L}],
MaxCellMeasure -> {1 -> h1}
];
Γ2 = DiscretizeRegion[
Line[{{-a, a/2}, {-a, a}, {a, a}, {a, -a}, {-a, -a}, {-a, -a/2}}],
MaxCellMeasure -> {1 -> h2}
];
Γ = RegionUnion[Γ1, Γ2];
(* triangle refinement function*)
cf = With[{h1 = h1, h2 = h2, a = a},
Compile[{{c, _Real, 2}, {area, _Real, 0}},
If[Max[Abs[c]] > 1.4 a, area > h1^2/2, area > h2^2/2.1],
CompilationTarget -> "C"
]
];
Ω = Module[{inst, outInst},
inst = TriangleCreate[];
TriangleSetPoints[inst, MeshCoordinates[Γ]];
TriangleSetSegments[inst,
MeshCells[Γ, 1, "Multicells" -> True][[1, 1]]];
outInst = TriangleTriangulate[ inst, "pq30aYY", "TriangleRefinement" -> cf];
MeshRegion[TriangleGetPoints[outInst], Triangle[TriangleGetElements[outInst]]]
];
numdof = MeshCellCount[Ω, 0];
pts = MeshCoordinates[Ω];
edges = MeshCells[Ω, 1, "Multicells" -> True][[1, 1]];
edgelookuptable = SparseArray[
Rule[
Join[edges, Transpose[Reverse[Transpose[edges]]]],
Join @@ ConstantArray[Range[Length[edges]], 2]
],
{numdof, numdof}
];
Creating a function to look up the points and edge of Γ1 and Γ2 in Ω.
LookupPoints[pts_, edgelookuptable_, Γpts_] :=
Module[{vertices, edgelist, numdof, nf},
numdof = Length[pts];
nf = Nearest[pts -> Automatic];
vertices = Flatten[nf[Γpts]];
edgelist = DeleteDuplicates@Sort@Extract[
edgelookuptable,
Times[
KroneckerProduct[#, #] &@
SparseArray[Partition[vertices, 1] -> 1, {numdof}],
edgelookuptable
]["NonzeroPositions"]];
{vertices, edgelist}
];
Looking up the points and edge of Γ1 and Γ2 in Ω.
{Γ1vertexlist, Γ1edgelist} = LookupPoints[pts, edgelookuptable, MeshCoordinates[Γ1]];
{Γ2vertexlist, Γ2edgelist} = LookupPoints[pts, edgelookuptable, MeshCoordinates[Γ2]];
HighlightMesh[Ω, {{Line[edges[[Γ1edgelist]]]}, {Line[edges[[Γ2edgelist]]]}}]
Using the low-level FEM-tools to set up the constrained linear system and solving it with LinearSolve.
(*setting up as much of our PDE with NDSolve`FEM` as possible*)
Needs["NDSolve`FEM`"]
Ωdiscr = ToElementMesh[Ω, "MeshOrder" -> 1, MaxCellMeasure -> ∞];
ClearAll[u];
vd = NDSolve`VariableData[{"DependentVariables", "Space"} -> {{u}, {x, y}}];
sd = NDSolve`SolutionData[{"Space"} -> {Ωdiscr}];
cdata = InitializePDECoefficients[vd, sd,
"DiffusionCoefficients" -> {{-ν IdentityMatrix[2]}},
"MassCoefficients" -> {{1}},
"LoadCoefficients" -> {{1}}
];
bcdata = InitializeBoundaryConditions[vd, sd, {{NeumannValue[0., True]}}];
mdata = InitializePDEMethodData[vd, sd];
dpde = DiscretizePDE[cdata, mdata, sd];
dbc = DiscretizeBoundaryConditions[bcdata, mdata, sd];
{load, stiffness, damping, mass} = dpde["All"];
DeployBoundaryConditions[{load, stiffness}, dbc];
(*stripping the Dirichlet vertices; being no boundary points,
Mathematica wouldn't do that for me through the FEM interface*)
plist = Complement[Range[numdof], Γ2vertexlist];
(*setting up the saddle point system*)
ξ = SparseArray[{Total[mass][[plist]]}];
A = ArrayFlatten[{
{stiffness[[plist, plist]], Transpose[ξ]},
{ξ, 0.}
}];
b = Join[Flatten[load[[plist]]], {1.}];
(*solving...*)
{u1, λ} = Through[{Most, Last}[LinearSolve[A, b]]]; // AbsoluteTiming
(*adding Dirichlet DOFs*)
u = ConstantArray[0., {numdof}];
u[[plist]] = u1;
ufun = ElementMeshInterpolation[{Ωdiscr}, u, InterpolationOrder -> 1];
Finally, the long awaited plots of the solution:
Plot3D[ufun[x, y], {x, y} ∈ Ωdiscr,
NormalsFunction -> None,
AxesLabel -> {"x", "y", "u"}
]
ContourPlot[ufun[x, y], {x, y} ∈ Ωdiscr, Contours -> 40]
Actually, things are not that easy.
The linear equations do not impose restrictions on the direction of matter transfer. If I am not completely mistaken, they would allow matter to disappear in free space and to be reinjected through $\varGamma_2$. In order to prohibit that, one has to add an outflow constraint of the form
$$N_2 u \leq 0 \quad \text{almost everywhere on $\varGamma_2$}.$$
Now this problem is a linear complementarity problem. While methods to tackle it in the finite dimensional discretization exist, I absolutely don't know whether it is well defined in the infinite-dimensional setting.
In general, $N_2 u$ can only be expected to be a member of $H^{-1/2}(\varGamma_1) := (H^{1/2}(\varGamma_1))'$ and one has to think about how $N_2 u \leq 0$ should be formulated. Maybe as
$N_2 u \in \operatorname{pol}(K)$
being a member of the polar cone of $K = \{v \in H^{1/2}(\varGamma_1)| v\geq 0\}$:
$$N_2 u \in \operatorname{pol}(K) := \{ \xi \in (H^{1/2}(\varGamma_1))' | \forall v \in K\colon \langle \xi, v \rangle \leq 0 \}. $$
(The cone $K$ is easily shown to be closed and convex.)
Looking at the solution of the example at hand, $N_2 u$ appears to be nonpositive on all of $\varGamma_2$. So this complicated setup might not be needed. But notice that $N_2 u \leq 0$ actually says only that the bulk outflow is dominant; single particles could still travel backwards. So, I am not entirely sure whether I modelled the one-way teleporter correctly.
