Im trying to solve a system of two PDE (im(t,x) and ia(t,x)) that are dependent on time and distance. I’m using my own anisotropic mesh (thanks to @TimLaska) that is aimed to represent an interface between a membrane and a liquid. The interface is located at L/2 where L = thickness of the system. “im” takes positive values in the interval 0 <= x <= L/2 and is equal to 0 at x > L/2. “ia” behaves similarly.
The problem that I’m encountering is that when I specify a Neumann boundary value at L/2, Mathematica doesn’t find this value in the mesh and the Neumann value is effectively ignored.
These are the functions for creating anisotropic meshes (kindly provided by @TimLaska in one of my previous questions)
(*Import required FEM package*)
Needs("NDSolve`FEM`");
(*Define Some Helper Functions For Structured Meshes*)
pointsToMesh(data_) := MeshRegion(Transpose({data}), Line@Table({i, i + 1}, {i,Length(data)- 1}));
unitMeshGrowth(n_, r_) := Table((r^(j/(-1 + n)) - 1.)/(r - 1.), {j, 0, n - 1})
meshGrowth(x0_, xf_, n_, r_) := (xf - x0) unitMeshGrowth(n, r) + x0
firstElmHeight(x0_, xf_, n_, r_) := Abs@First@Differences@meshGrowth(x0, xf, n, r)
lastElmHeight(x0_, xf_, n_, r_) := Abs@Last@Differences@meshGrowth(x0, xf, n, r)
findGrowthRate(x0_, xf_, n_, fElm_) := Quiet@Abs@ FindRoot(firstElmHeight(x0, xf, n, r) - fElm, {r, 1.0001, 100000},Method -> "Brent")((1, 2))
meshGrowthByElm(x0_, xf_, n_, fElm_) := N@Sort@Chop@meshGrowth(x0, xf, n, findGrowthRate(x0, xf, n, fElm))
meshGrowthByElm0(len_, n_, fElm_) := meshGrowthByElm(0, len, n, fElm)
flipSegment(l_) := (#1 - #2) & @@ {First(#), #} &@Reverse(l);
leftSegmentGrowth(len_, n_, fElm_) := meshGrowthByElm0(len, n, fElm)
rightSegmentGrowth(len_, n_, fElm_) := Module({seg}, seg = leftSegmentGrowth(len, n, fElm);
flipSegment(seg))
extendMesh(mesh_, newmesh_) := Union(mesh, Max@mesh + newmesh)
Here I define a couple of constants and create my own mesh
(* Constants *)
L = 0.1; dim = dia = 1.*^-6; kf = kr = 1; kr =
1./kf;
(* Mesh generation *)
seg1 = rightSegmentGrowth(L/2., 100, L/1000);
seg2 = leftSegmentGrowth(L/2., 100, L/1000);
totalseg = extendMesh(seg1, seg2);
rh = pointsToMesh@totalseg ;
crd = MeshCoordinates(rh);
mesh = ToElementMesh(crd);
And here is the simple code I use to find the solution alongside the error I get
(*PDE system*)
eqsim = {D(im(t, x), t) - dim D(im(t, x), x, x) == NeumannValue(kf*ia(t, x), x == L/2.) + NeumannValue(0, x == L/2.),im(0, x) == 0.1};
eqsia = {D(ia(t, x), t) - dia D(ia(t, x), x, x) == NeumannValue(kr*im(t, x), x == L/2.) + NeumannValue(0, x == L/2.), ia(0, x) == 0.1};
sol = NDSolve({eqsim, eqsia}, {im, ia},x (Element) mesh, {t, 0, 60}, Method ->{"MethodOfLines", "SpatialDiscretization" -> {"FiniteElement"}})
(*NDSolve::bcnop: No places were found on the boundary where x==0.05 was True, so NeumannValue(ia,x==0.05) will effectively be ignored.*)
(*NDSolve::bcnop: No places were found on the boundary where x==0.05 was True, so NeumannValue(1. im,x==0.05) will effectively be ignored.*)
(* NDSolve::bcnop: No places were found on the boundary where x==0.05 was True, so NeumannValue(0,x==0.05) will effectively be ignored.*)
(*General::stop: Further output of NDSolve::bcnop will be suppressed during this calculation.*)
I first thought that the discretisation in the mesh didn’t include the point L/2 but after double checking I found out that it’s there
Position(crd, L/2.)
(*{{100, 1}}*)
Similar problems to this have been reported when using Dirichlet conditions here DiscretizeRegion does not include the boundary specified in ImplicitRegion (10.1)
, here PeriodicBoundaryConditions: missing points (a simpler example)
and here Solving Laplace’s equation in 2D using region primitives
, but after reading the solutions provided I haven’t been able to find a solution to my issue.
Any help, feedback or advice is more than welcome (I’m using v 12.2.0.0).