I have a problem with solution of the 2D PDE system. It appears that solution does not match boundary condition at all.
I've tried using separately generated mesh, and it had not helped, just making error oscillation frequency larger, but not lowering the error amplitude.
Adding "Method -> {"ExplicitRungeKutta", "SpatialDiscretization" -> {"FiniteElement", "MeshOptions" -> {"MaxCellMeasure" -> 0.00001}}}" does not make any difference.
The region is a rectangle (I tried oval too) with two round holes. There is no such problem in the case of just one hole.
Do you have any ideas making it better?
system = {nxU[x, y] == D[HU[x, y], x] + (D[nyU[x, y], x, y] + D[nxU[x, y], x, x]),
nyU[x, y] == D[HU[x, y], y] + (D[nyU[x, y], y, y] + D[nxU[x, y], x, y]),
D[nyU[x, y], y] + D[nxU[x, y], x] == D[HU[x, y], y, y] + D[HU[x, y],x,x]};
FuncLst = {nxU[x, y], nyU[x, y], HU[x, y]};
n0=0.3;hIn0=-0.5;rOut=30;r0=1.;Centre1={-10,0};Centre2={10,0};
RegionOut=Rectangle[{-rOut,-rOut/2},{rOut,rOut/2}];
BoundaryOut=RegionBoundary[RegionOut];
RegionIn1=Disk[Centre1,r0];
RegionIn2=Disk[Centre2,r0];
RegionIn=RegionUnion[RegionIn1,RegionIn2];
BoundaryIn1=RegionBoundary[RegionIn1];
BoundaryIn2=RegionBoundary[RegionIn2];
BoundaryIn=RegionBoundary[RegionIn];
RegionAll=RegionUnion[RegionDifference[RegionOut,RegionIn],BoundaryIn];
n0func1x[x_,y_]=n0 (x-Centre1[[1]])/r0;
n0func1y[x_,y_]=n0 (y-Centre1[[2]])/r0;
n0func2x[x_,y_]=n0 (x-Centre2[[1]])/r0;
n0func2y[x_,y_]=n0 (y-Centre2[[2]])/r0;
BndCond={DirichletCondition[nxU[x,y]==0,{x,y}\[Element]BoundaryOut],
DirichletCondition[nyU[x,y]==0,{x,y}\[Element]BoundaryOut],
DirichletCondition[HU[x,y]==0,{x,y}\[Element]BoundaryOut],
DirichletCondition[nxU[x,y]==n0func1x[x,y],{x,y}\[Element]BoundaryIn1],
DirichletCondition[nyU[x,y]==n0func1y[x,y],{x,y}\[Element]BoundaryIn1],
DirichletCondition[nxU[x,y]==n0func2x[x,y],{x,y}\[Element]BoundaryIn2],
DirichletCondition[nyU[x,y]==n0func2y[x,y],{x,y}\[Element]BoundaryIn2],
DirichletCondition[HU[x,y]==hIn,{x,y}\[Element]BoundaryIn1],
DirichletCondition[HU[x,y]==hIn,{x,y}\[Element]BoundaryIn2]};
solPDE = NDSolve[Join[system, (BndCond /. hIn -> hIn0)], Join[FuncLst],{x, y} \[Element] RegionAll][[1]];
nxUs[x_, y_] = Evaluate[nxU[x, y] /. solPDE];
xs[t_] = Centre1[[1]] + Cos[t];
ys[t_] = Centre1[[2]] + Sin[t];
Plot[{nxUs[xs[t], ys[t]], n0func1x[xs[t], ys[t]]}, {t, 0, 2 \[Pi]}]
FunLst
when I think you meanFuncLst
. $\endgroup$ – Michael Seifert May 24 '18 at 20:46